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Abstract 

High-resolution images of the solar surface show a granulation pattern of hot rising and cooler downward-sinking 
material - the top of the deep-reaching solar convection zone. Convection plays a role for the thermal structure of the 
solar interior and the dynamo acting there, for the stratification of the photosphere, where most of the visible light is 
emitted, as well as for the energy budget of the spectacular processes in the chromosphere and corona. Convective 
stellar atmospheres can be modeled by numerically solving the coupled equations of (magneto)hydrodynamics and 
non-local radiation transport in the presence of a gravity field. The C05BOLD code described in this article is 
designed for so-called "realistic" simulations that take into account the detailed microphysics under the conditions 
in solar or stellar surface layers (equation-of-state and optical properties of the matter). These simulations indeed 
deserve the label "realistic" because they reproduce the various observables very well - with only minor differences 
between different implementations. The agreement with observations has improved over time and the simulations are 
now well-established and have been performed for a number of stars. Still, severe challenges are encountered when 
it comes to extending these simulations to include ideally the entire star or substellar object: the strong stratification 
leads to completely different conditions in the interior, the photosphere, and the corona. Simulations have to cover 
spatial scales from the sub-granular level to the stellar diameter and time scales from photospheric wave travel times 
to stellar rotation or dynamo cycle periods. Various non-equilibrium processes have to be taken into account. Last 
but not least, realistic simulations are based on detailed microphysics and depend on the quality of the input data, 
which can be the actual accuracy limiter. This article provides an overview of the physical problem and the numerical 
solution and the capabilities of C05BOLD, illustrated with a number of applications. 
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1. Introduction 

In the core of the Sun, fusion of hydrogen to helium releases energy which is transported outward, first by radiation 
only, but further out primarily by convection in the outer 30 % of the radial distance to the solar surface. Most of this 
energy is emitted in the form of radiation in the photosphere which is the bottom layer of the solar atmosphere. 
Furthermore, a small part of the energy is carried by waves and by magnetic fields, powering the dramatic phenomena 
visible in the solar chromosphere and corona. In more massive and further evolved stars, the internal structure is more 
complex, with several shells where nuclear burning takes place and multiple convection zones. 

The relatively thin solar photosphere (about 0.1 % of the solar radius) therefore plays an important role for the 
inner as well as for the outer layers of the Sun. The analysis of solar and stellar spectra can reveal surface properties 
and the chemical composition, and allows us to draw conclusions about the internal structure and evolutionary status. 
For this purpose, physical models of stellar atmospheres with a realistic treatment of both radiation and convection 
are essential. 
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The classical analysis relies on one-dimensional (ID) stationary model atmospheres (in most cases only the pho- 
tosphere plus the very top layers of the surface convection zone), where the average convective energy flux F conv is 
computed from the so-called local mixing-length theory []]], J^, a heuristic recipe which assumes that Fconv can 
be determined from local properties of the stratification. In the framework of this "theory", the mean thermal structure 
of a convective stellar atmosphere is found by the requirement that the sum of radiative and convective flux equals the 
total stellar flux, F rac j + F CO!W = cr r^ ff , at all depths. State-of-the-art radiative-convective equilibrium models of solar 
and stellar atmospheres have been constructed with the classical model atmosphere codes ATLAS |H|^], MARCS 
|S0]j an d PHOENIX flm, to name the most prominent examples. 

A severe drawback of these models is that the efficiency of the convective energy transport is controlled by a free 
parameter, the mixing-length parameter (Tmlt, which is of the order unity but a priori unknown. Therefore, amr 
must be calibrated against observations. Unfortunately, different observables require different values of ojmlt 11011 . 
The best fit of the Balmer line profiles of solar-type stars is achieved with (Xmlt ~ 0.5 111 ill , while continuum colors 
are better reproduced with ckmlt in the range 1-2, depending on the considered wavelength range fll2[]. The standard 
stellar-evolution calibration based on matching the current solar parameters calls for q-mlt ~ 2 lll3ll . ffl . This 
disparity indicates that the underlying theoretical description is inadequate. In fact, the solar photosphere is neither 
homogenous nor static, since it is influenced by the very top of the convection zone and shows a granular pattern of 
bright upflow regions surrounded by darker intergranular lanes of downflowing material, with a spatial scale of about 
1 Mm (10 6 m) and evolving on a time scale of minutes (see Fig.[T]for snapshots from two C05BOLD simulations of 
the solar granulation). This motivated various efforts to overcome the limitations of the ID classical atmospheres and 
to develop instead self-consistent, parameter-free hydrodynamical models of stellar surface convection, accounting 
for the fact that convection is a non-local, time-dependent, and intrinsically three-dimensional phenomenon. 

Early idealized numerical simulations of convection under stellar-like conditions had to resort to severe simplifi- 
cations (stationary 2D solutions on coarse grids) and could only deliver qualitative results: Latour et al. fl5i[l6ll and 



and Nelson 



Toomre et al. 117 [used anelastic modal equations to study surface convection in A-type stars. Musman and Nelson 



1911 investigated convection in the Sun and some other stars with a similar method. Chan and Wolff 



12011 developed a code based on the alternating direction implicit (ADI) method for the calculation of compressible 
convection. Hurlburt et al. 12111 carried out simulations of compressible solar convection extending over multiple scale 
heights. Steffen et al. [22] took (non-local) radiation transfer into account in their 2D simulations of compressible solar 
convection. 

The first realistic simulations of solar granulation were performed by Nordlund 12311 and included three-dimensional 
(3D) time-dependent hydrodynamics (but anelastic and with moderate spatial resolution) and non-local radiative en- 
ergy transfer, already then with a simple treatment of the frequency-dependence of the opacities. Hand-in-hand came 
the a posteriori detailed spectrum synthesi s by Dravins et al. 12411 . Other 3D convection simulations relinquished the 
treatment of radiation transfer (25, 26[ 27, 28, 29], Current radiation hydrodynamic codes of various groups use sim- 
ilar basic techniques - in a significantly refined way (compressible hydrodynamics, more grid points, more opacity 
bins, larger computational domains, magnetic fields, a chemical reaction network, dust, etc.). For example, Stein and 
Nordlund have carried out radiation hydrodynamics (RHD) simulations with 2016 x 2016 x 500 grid points with a 
spatial resolution of 24 km in the horizontal direction and 12-80 km in the vertical direction. Asplund et al. |3(i 31] 
have computed chemical abundances using high spatial resolution and accurate radiative transfer. 

Just like classical stellar atmospheres, the non-magnetic hydrodynamical models are characterized by the average 
total energy flux per unit area and time (effective temperature, r e ff), surface gravity g, and chemical composition. 
But, in contrast to the mixing-length models, there is no longer any free parameter to adjust the efficiency of the 
convective energy transport. Similarly, the fudge parameters micro- and macroturbulence, that have to be introduced 
in ID model atmospheres to match synthetic and observed shapes of spectral lines, are replaced by the self-consistent 
hydrodynamical velocity field of the 3D simulations. However, one has to keep in mind that the simulations are 
characterized by a large number of numerical parameters, e.g., the spatial resolution of the numerical grid, the size 
of the computational domain, the formulation of boundary conditions, and the parameters related to the numerical 
schemes for solving the hydrodynamical and radiation transport equations. Of course, the hope is that the simulation 
results become essentially independent of the choice of these numerical parameters, once a sufficiently high spatial, 
angular, and frequency resolution is achieved. 

Hydrodynamical model atmospheres are not only computed for the Sun but also for other stars, and are comple- 
menting and increasingly replacing classical ID atmosphere models. Important applications of convection simulations 
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Figure 1 : Emergent continuum intensity at a wavelength of A = 500 nm, synthesized for a snapshot from a high-resolution (400 X 400 X 300) 
C05BOLD RHD simulation (left, see Sect.gT} and for a (286 X 286 X 266) C05BOLD RMHD model (right, cf. Fig.[l4j, each representing a small 
patch of solar surface granulation. The imposed average magnetic flux of the RMHD model is (B z ) = 50 G. The arrows represent streamlines that 
follow the horizontal velocity on the surface of optical depth unity, i.e., at the bottom of the solar photosphere. 



with C05BOLD and its predecessor include the accurate spectroscopic determination of solar and stellar chemical 



abundances and isotopic ratios [e.g., |32|,|33l|34|], the theoretical calibration of the mixing-length parameter 13511 . the 



study of convective overshoot and mixing processes in stellar envelopes H36H . and the excitation of waves by turbulent 
convective flows 137,38, 39| . 

The presence of magnetic fields results in a wide range of additional complex 3D phenomena. Small-scale concen- 
trations of magnetic flux lead to enhanced radiative losses, both in the photosphere and in the chromosphere. On the 
other hand, large-scale magnetic structures can inhibit the convective energy flux and produce the well-known dark 
sunspots. The interaction of convection and magnetic fields can be modeled in the framework of (ideal) magneto- 
hydrodynamics (MHD). 

In the purely hydrodynamical simulations described above, the resulting mean flow is determined only by the 
prescribed physical quantities T e g , g, and the assumed chemical composition, and is largely independent of the for- 
mulation of the boundary conditions and details of the initial configuration. This is no longer true for the more 
complex simulations of solar magnetoconvection. In this case, the presence of a magnetic field implies more freedom 
in setting up the problem: the initial configuration of the magnetic field and the magnetic boundary conditions have 
to be designed for the particular problem under consideration. In many studies, the magnetic field is assumed to be 
vertical at the upper and lower boundaries, such that the horizontally averaged magnetic flux is fixed at a prescribed 
value. For example, {B : ) — 50 G for the C05BOLD MHD simulation shown in Fig.[T] which is representative of the 
least magnetic solar-surface areas, the so-called quiet-Sun internetwork regions. The velocity arrows in this figure 
show that the flow converges towards the dark intergranular lanes, where cool gas returns into the solar convection 
zone. This flow also leads to a concentration of magnetic flux in the downflow lanes, where is is visible as bright 
knots or elongated features (e.g., near x = 5.5 Mm, y = 1.5 Mm in Fig.Q] right panel). 

Early 2D MHD simulations of solar convection, which include radiative transfer were presented by Grossmann- 
Doerth et al. [40] based on a adaptive moving finite element code, by Steiner et al. 114 111 with a finite-volume code 
based on automatic adaptive mesh refinement for MHD, described in Steiner et al. B42I1 . and by Atroshchenko and 
Sheminova [43] who used a method of approximate Eddington factors for the radiative transfer. 

To our knowledge, the first realistic three-dimensional radiation magnetohydrodynamic (RMHD) simulation of 
stellar magnetoconvection was presented by Nordlund 14411 . Nordlund et al. 14511 give a review on solar surface 
convection including results on magnetoconvection. Early two-dimensional MHD simulations of stellar magneto- 
convection, which dispense with detailed radiative transfer include Galloway and Weiss fi^l . Deinzer et al. 1471. 
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Hurlburt and Toomre (1988) [48], Weiss et. al. (1990) |49J], and Fox et al. (1991) |5C 

The pioneering work of Nordlund and collaborators was only recently followed up by others, also working in 
three spatial dimensions. Examples include Hansteen and Gudiksen |51f and Gudiksen et al. 05211 with the Bifrost 
code, Schaffenberger et al. J53_ with C05BOLdH] and Vogler et al. ll54ll with the MURaM codeJland more recently, 
by Heinemann et al. 1 55 ] using the Pencil code0Jacoutot et al. |[56ll with a code named SolarBox, developed by A. 
Wray, and Muthsam et al. 15711 with the Antares code. Recent impressive large-scale 3D RMHD simulations include 
the supergranulation-size magnetoconvection simulations by Stein et al. [58], using a variant of the STAGGER code 
of Nordlund and Galsgaard |59], the simulations of sunspots and solar active regions described in Cheung et al. IrSoTl 
and in Rempel et al. 116 ill , both works using the MURaM code, as well as the exploratory MHD models that span 
the entire solar atmosphere from the upper convection zone to the lower corona by Hansteen et al. [62], |63_, and 
Martmez-Sykora et al. „64Tl . based on Bifrost or an extended version of the STAGGER code. 

Other three-dimensional simulations of stellar magnetoconvection use approximations to the radiation transfer, 
like Abbett J65ll . Abbett and Fisher „66ll . and Isobe et al. |67[. Important results of solar magnetoconvection in three 
spatial dimensions were also obtained bysimply replacing the radiation transfer with heat conduction, e.g., by Weiss 
et al. l68ll . Tobias et al. ll69ll . Cattaneo I70h . Ossendrijver et al. |71], or Cattaneo et al. |72[. For other applications, 
radiative exchange or heat conduction is not as critical as for convection, e.g., for the rise of buoyant magnetic flux 
tubes. Such simulations were carried out, e.g., by Archontis et al. l73ll with the STAGGER code J59ll or by Cheung 
et al. Q with the Flash codeQ 

Simulations of global stellar convective dynamos have been started by Glatzmaier IT75ll . More recent global MHD 
simulations of stellar convection include Browning et al. lf76ll with the ASH-code |77] and Dobler et al. 17811 with the 
Pencil code. Ziegler |79] applied the Nirvana codqjto the problem of core collapse and fragmentation of a magnetized 
protostellar cloud. 

Further MHD codes for potential application to realistic stellar convection simulations, which have been developed 
in an astrophysical context are the A-MAZE code@ the Enzo codeQ the VAC codeJI or the Zeus code@ for a non 
exhaustive list. 

For reviews on solar magnetoconvection see Nordlund et al. fi^ll . Nordlund and Stein l80ll . Carlsson jilll . Steiner 



2. Basics 

2.1. Basic considerations about convective scales 

Ideally, hydrodynamical models of stellar convection should comprise the entire convection zone in a spherical 
shell with sufficient spatial resolution, and should cover all relevant time scales. In general, such a global approach is 
not feasible, however, for the reasons outlined in the following basic considerations. 

2.1.1. Spatial scales 

Presently, realistic models of stellar convection are restricted to a small representative volume located near the 
surface, including both the top layers of the convection zone and the photosphere, where most of the stellar radiation 
is emitted. In this context, it is important to realize that convection is driven by entropy fluctuations generated near 
the surface by radiative cooling. The deeper layers approach an adiabatic mean state and have little direct influence 
on the small-scale granular flows at the surface. For this reason, it is possible to obtain physically consistent ab initio 
models of stellar surface convection from local-box simulations that cover only a small fraction of the geometrical 



1 See http://www.astro.uu.se/~bf/co5bold_main.html. http://www.co5bold.com 
2 See http://www.mps.mpg.de/projects/solar-mhd/muram_site/code.html 
3 See http://www.nordita.org/software/pencil-code/ 
4 See http://flash.uchicago.edu/website/home/ 
5 See http://nirvana-code.aip.de/ 

6 See http://www.the-a-maze.net/people/folini/research/a_maze/a_maze.html 
7 See http://lca.ucsd.edu/portal/software/enzo 
8 See http://grid.engin.umich.edu/~gtoth/VAC/ 
9 See http://lca.ucsd.edu/portaI/codes/zeusmp2 
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Figure 2: Reynolds number Re (left) according to Eq. f^, and Prandtl number Pr (right) computed from Eqs. (4} and UOt with k = ko = lOn/Hp, 
as function of radius in the envelope and atmosphere of the Sun, using the solar model by Christensen-Dalsgaard et al. 1131 . In addition, the upper 
(dashed) curve in the right panel refers to the Mach number, Ma = v c /c s . The vertical dashed line marks the bottom of the solar convection zone. 



depth of the whole convection zone. Since the lower boundary is thus located right inside the convection zone where 
the total stellar luminosity is entirely carried by the convective flow, it is essential to employ an open lower boundary 
condition that impedes the flow as little as possible (details are given in Sect. l3.2.TT i. 

As a typical example, let us consider a local -box simulation of the solar granulation measuring L x x L y — 10 Mm 
xlO Mm in the horizontal directions with periodic lateral boundary conditions in x and y. In the vertical direction, 
open boundaries are imposed, and the extension of the box is assumed to be L : = 4 Mm, with L~ w 3 Mm (A In P « 7) 
below and L* « 1 Mm (AlnP w 8) above the optical surface, where A In P is the number of gas pressure e-foldings. A 
box of this size covers only 1.5% of the total depth of the solar convection zone, but is large enough to accommodate 
several surface convection cells called granules (cf. Fig.[T]i, ensuring that the periodic boundary conditions do not have 
a critical influence on the resulting flow pattern. The minimum spatial resolution of the numerical grid is set by the 
requirement to cover one pressure scale height by at least 10 grid cells. In the following, we assume that a typical 
grid comprises N x X N y x N z = 250 x 250 x 200 cells, where the horizontal cell size is constant (Ax = Ay * 40 km), 
while the vertical cell size increases with depth (in proportion to the local pressure scale height H p , see below) from 
about 10 km at the surface to about 50 km near the bottom of the computational domain (for some actual examples 
see Table [TJ. 

It is well known that the convective envelope of the Sun is characterized by very large flow Reynolds numbers, Re. 
Based on the standard solar model of Christensen-Dalsgaard et al. lfl3ll . we have evaluated this dimensionless number 
locally as 

Re = , (1) 

v 

where H p = -(dlnP/dz) -1 is the local pressure scale height (e-folding length of the gas pressure P), v c is the charac- 
teristic convective velocity according to classical mixing-length theory yj,|2|], and v is the microscopic (atomic plu s 
radiative) kinematic viscosity, v = {r\ a + rj r )/p, with rj a and rj r calculated according to Spitzer 18311 and Thomas 1184m 1 . 
respectively. The depth dependence of Re in the solar envelope is displayed in the left panel of Fig. [2] showing that 
Re > 10 10 in the entire convection zone. This implies that the flow is highly turbulent wherever convection occurs 
(see however 18511 ). The turbulent kinetic energy is dissipated into heat at the Komolgorov microscale, £ « H p Re~ 3 ^ 4 , 
which varies between 0.05 and 10 cm from the top to the base of the solar convection zone. Clearly, the spatial resolu- 
tion of the numerical simulations sketched above is insufficient by more than 6 orders of magnitude to properly resolve 
the complete turbulent cascade. All realistic stellar convection simulations therefore follow the so-called large-eddy 
approach, where only the largest flow structures, including the driving scales, are resolved, and the small-scale kinetic 
energy is dissipated at the grid scale, either by the numerical scheme or by a subgrid-scale model. Consequently, the 
effective numerical viscosity in such models is at least 8 orders of magnitude larger than in reality. 

In addition to the Reynolds number, the properties of the flow are further characterized by the (dimensionless) 
Prandtl number, 

Pr = - , (2) 
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Figure 3: Left: Kelvin-Helmholtz time scale tkh (upper set of curves, computed from Eqs. (5), (6)), convective turnover time scale Turnover (middle 
set of curves, computed according to Eqs. 0, {8)), and the CFL time scale tcfl (lower curve, Eq. j9j), as a function of radius in the solar envelope, 
based on the solar model by Christensen-Dalsgaard et al. fiUl . The vertical dashed line in the left panel indicates the bottom of the solar convection 
zone. The right panel zooms into the upper 5% of the convection zone, and shows in addition the radiative time scale T ra( j calculated from Eq. IIP) 
with k = k a = 10 it /H p (dotted). 

the ratio of the coefficients describing the diffusion of momentum, v, and heat, x- In the stellar interior and atmosphere, 
heat transfer is dominated by radiation, which in the optically thick layers can be described as a diffusion process. 
The radiative diffusivity is given by 



(<x: Stefan-Boltzmann constant, T: temperature, k: radiative opacity per unit mass, p: mass density, c v : specific 
heat at constant volume). Pr depends only on the thermodynamic state of the stellar gas. In the optically thin layers 
(photosphere), radiative heat exchange cannot be described as a diffusion process, and hence the definition of Pr via 
Eqs. (f2| and (O is no longer meaningful. Instead, Pr can be defined more generally as 



the ratio of radiative time scale (f ra( j, see Eq. (TToT> below) to viscous time scale (fv. = vk 2 ). However, the Prandtl 
number then becomes a function of wavenumber k for optically thin conditions. In the solar convection zone and 
atmosphere, Pr ranges between 10~ 4 and 10~ 10 (see Fig.|2j right panel), indicating that the radiative energy diffusion 
is much more efficient than the viscous diffusion of momentum, in other words, the dynamical lifetime of a turbulent 
vortex is much longer than its thermal relaxation time. 

In large-eddy simulations, diffusion is provided by an explicit artificial viscosity and/or by the numerical advection 
scheme, which leads to a diffusive cutoff at the scale of the grid resolution. In general, the effective viscosity depends 
on the grid resolution Ax, and on the wavenumber (and amplitude) of the local velocity perturbation. For small-scale 
structures close to the grid resolution, the coefficients characterizing the numerical diffusion of momentum, v, and 
heat, x, are of similar size, and hence the Prandtl number is of the order unity, Pr = v/x ~ 1, as long as the radiative 
diffusivity is much smaller than the numerical one, % <K x ~ "c Ax. This condition always holds in the bulk of the 
solar convection zone (assuming Ax w H p /1Q). On the other hand, the effective artificial/numerical diffusion can 
be significantly smaller for well resolved smooth structures, such that x > X > an d Pr ~ v/x < 1 • This is especially 
true for the near-surface layers where the radiative diffusivity is high. Large-eddy simulations of solar-type surface 
convection can therefore achieve moderately low Prandtl numbers, in the sense that the physical radiative energy 
transport dominates over numerical diffusion of heat. In the bulk of the convection zone, however, the radiative 
diffusivity is too low, and hence numerically Pr * 1 on all resolved scales. 

2.1.2. Time scales 

The time span covered by a numerical convection simulation must be sufficiently long to ensure that the whole 
structure contained in the computational box can reach a thermally relaxed state. Thermal relaxation by radiative 
diffusion proceeds on the Kelvin-Helmholtz time scale, defined as the thermal energy content (per unit area) divided 




(4) 



by the total energy flux: tkh = E t h/F to t- In Fig. [3] we show the depth-dependence of tkh in the solar convection zone 
(upper curves), computed as 



r {1) ( 



4tt r 2 pc„T H„ 

. (5) 



and 



T ™W= T- f° pc P Tr 2 dr , (6) 

^0 Jr 



respectively (R Q and L denote the solar radius and luminosity). Both expressions give essentially identical results, 
indicating that at a depth of z = -3 Mm below the solar surface, tkh ~ 10 6 s or « 280 h. Fortunately, it turns out 
that relaxation is significantly faster than expected from this estimate. Since the energy flux is carried by convection, 
a few convective turnover times are sufficient to establish a self-consistent equilibrium state. The convective turnover 
time scales, calculated as 

t {1) (r) = ^ (7) 

turnover V / » v / 



and 



Jr 



^rnoverW = | f & , (8) 



respectively, are also shown in Fig.[3](left, middle curves). The plot shows that T turnover ~ 2000 s at z = -3 Mm, a factor 
500 smaller than tkh- Roughly, the simulation needs to be advanced for about 10 turnover times, f s j m » 10 r turnover w 
20 000 s to obtain a relaxed model. This number has to be related to the numerical time step At applicable to the 
hydrodynamics scheme. The well-known Courant-Friedrichs-Lewy (CFL) condition for the stability of an explicit 
numerical method states that At < tcfl, where tcfl is given by the travel time of the fastest wave across a grid cell. 
For the present non-magnetic example we can use the approximation 

Ax H p 

tcflO) = ~ 7777 7 , (9) 

(c s + v c ) 10 (c s + v c ) 

where c s is the adiabatic sound speed. Evaluation of Eq. (0 shows that tcfl ~ 1 s in the upper layers (see Fig.0. 

The numerical time step is not only limited by the CFL condition. In addition, Af must be smaller than the 
characteristic radiative time scale T ra( j that rules the decay of local temperature perturbations at the smallest possible 
spatial scale (wavenumber ko = 10n/H p ). To a good approximation, T ra d can be calculated as 

T ^ (r) = _£L_( 1 + 3 ^ = M-^ + r 2 ) . do) 



16ctkT 3 \ k 2 ) x \3p 2 K 2 

which is valid in both optically thick and thin regions |86[ 87]. As illustrated in Fig. [3] (right), r m d(ko) reaches a 



sharp local minimum of » 0.2 s close to the optical surface. The time step of the numerical simulation is thus set by 
the radiative time scale, Af < 0.2 s, and the total number of required time steps is N t = f s im/Af ~ 10 5 . Assuming 
for reference a processor that can update N c = 10 6 grid cells per CPU second, the total CPU time required for this 
standard simulation would be (A^ x N y x N z x N,)/N c ~ 10 12 /10 6 s or about 12 days, which is well feasible even 
without a high degree of parallelization. However, it is also clear that much larger models (e.g., 10 times better spatial 
resolution in each direction) are out of reach without massive parallelization. 

As an example, consider the solar supergranulation which has a typical horizontal scale of 20-30 Mm. Numerical 
simulations of this phenomenon thus require a horizontal cross section of at least 100 x 100 Mm 2 . Since the spatial 
resolution cannot be reduced much if the granular scale still is to be resolved, such a horizontally extended simulation 
would take a factor 100 more CPU time than the standard case outlined above. In addition, the simulation box would 
need to be extended to deeper layers for this kind of modeling. Assume that the lower boundary is moved from a depth 
of z = -3 Mm to z = -20 Mm, which means extending the model by about 6 more pressure scale heights. Adding 
100 grid cells in the vertical direction could be sufficient to cover the extra 17 Mm. In terms of computing time, these 
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extra cells are relatively cheap, because radiative transfer can be treated by the diffusion approximation in these deep 
layers. Note, however, that keeping the horizontal resolution at Ax = Ay ss 40 km to resolve the granulation at the 
surface, the aspect ratio of the cells near the bottom of the box becomes rather extreme, Ax/Az ~ 1/10. But the real 
problem is that the turnover time increases by a factor 20. Since At is set by the surface layers, the number of time 
steps increases by the same factor. In summary, a supergranulation simulation will take roughly a factor 2000 more 
time than a standard granulation model, about 65 years of CPU time. With massively parallel computers, such models 



are becoming marginally feasible (cf. 15811 ) 



2.1.3. Global convection simulations 

Simulations of the entire solar convection zone are much more expensive: the turnover time increases by another 
factor 100, while the surface area is about 600 times larger with respect to the above supergranulation model. In terms 
of the numbers quoted above, such a global convection simulation, which ideally should be carried out in a rotating 
spherical coordinate system, would take of the order of 4 million years of CPU time, but still would cover only one 
year of solar time. In order to study the solar magnetic dynamo action, it would certainly be desirable to run the 
simulation over several 22-year cycles, say a period 100 solar years, which is equivalent to 400 million CPU years. 

Since the surface layers set the numerical time step and spatial resolution, the computational cost can be much 
reduced by restricting the simulations to the deeper layers of the convection zone: here the flow Mach number is small 
(see Fig.|2]i, and the so-called anelastic approximation can be employed to avoid the time step limitation by the CFL 
condition; moreover, the radiative time step is very large (see Fig. [3} and does not impose any additional limitation. 
This approach has been adopted in the global simulations of the solar convection zone with the ASH-code by Brun 
et al. 08811 . However, the direct link between model and observation is necessarily broken in such kind of modeling. 

While realistic simulations of global solar convection remain phantasmal, prospects can be better for other type of 
stars: realistic global star-in-a-box simulations have already been performed successfully for red supergiants, where 
only a few huge convection cells occupy the surface of the star (see Sect. l4.7l i. 

2.1.4. From the upper convection zone to the lower corona 

The essential physics necessary for realistic simulations of solar surface convection includes compressible hy- 
drodynamics describing transonic flows of a partially ionized gas in a gravitationally stratified atmosphere, coupled 
with non-local, frequency-dependent radiative energy exchange. In the subsurface layers, the flow becomes strongly 
subsonic and can be described in the anelastic approximation, while the radiative transfer becomes local and can be 
treated by the gray diffusion approximation. In contrast, physics becomes more complicated when considering the 
outer solar atmosphere. 

Simulations comprising the chromosphere and lower corona must include magnetic fields. Since the magnetic field 
tends to form localized flux concentrations in the intergranular lanes (cf. Fig.Q] right panel), the spatial resolution of 
MHD simulations needs to be better than that of non-magnetic granulation models. In addition, the time step is 
dictated by the Alfven speed 

B 

VA=^= , (11) 



which can become much larger than the sound speed in places where the plasma-/? is low, i.e., where the magnetic 
field B is large and the density p is small. Typically, A?mhd ~ A?hd/100. 

The low density of the outer atmosphere has also consequences for the radiation transport. Since the collision 
frequency is reduced, the simplifying assumption of local thermodynamic equilibrium (LTE) tends to break down, 
and photon scattering becomes important. This implies that the source function is no longer a function of the local 
temperature, but depends also on the angle-averaged radiation field. In contrast to the photospheric absorption line 
spectrum, the chromospheric spectrum contains strong emission lines, which dominate the energetics in the chromo- 
sphere. Under these circumstances, the solution of the radiation transfer problem becomes very time consuming. 

Heat transfer by thermal conduction becomes important above gas temperatures of a few 10 4 K, i.e., in the transi- 
tion region and in the corona above [see, e.g.,|89l|93l. Thermal conduction is usually modeled by means of the Spitzer 
formula but can result in a significant increase of the computational costs. 

Further complications arise due to the fact that the ionization of hydrogen (and other elements) is no longer in 
thermal equilibrium in the low density regions, and cannot be obtained from precomputed look-up tables. Rather, 
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the degree of hydrogen ionization, and hence the electron density, has to be derived from the solution of the time- 
dependent rate equations of a multi-level atom, which poses severe challenges. 

For further discussion of these problems see Sect. 14. 5. Tl as well as Hansteen et al. JH^l, l63ll . Martmez-Sykora 
et al. and Gudiksen et al. iBll . 



2.2. Equations 

The hydrodynamics equations are expressed as conservation relations plus source terms for 
p, pu u pv 2 ,pv 3 , e, ot , 



(12) 



the mass density, the three momentum densities, and the total energy density (per volume), respectively. The coordi- 
nate axes are simply numbered, in this case and in the code itself. In some sections, we use the more standard notation 
x, y, and z, though. 

The three-dimensional hydrodynamics equations, including source terms due to gravity, are the mass conserva- 
tion equation 



dp d pvi d pv 2 d p v 3 

-f - + —7T + — f + —7T = , 

Ot OX\ 0X2 OX 3 

the momentum equation 
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and the energy equation 

dpe tot d(pe m +P)vi d(pe tot +P)u 2 8 (pe m +P) v 3 dF lr . dd 3F 2rad dF 3md 

+ 7. + 7. + 7. + -T. + — 7 + — " = . 



dt 



dx\ 



8x2 



dx 3 



dx\ 



dxo 



dx 3 



(13) 



(14) 



(15) 



Here F\ m A, ^2rad> ^3rad are the components of the radiative energy flux (see below). The gas pressure P is computed 
from the density p and the internal energy, e mt , via an equation of state, usually available to the program in tabulated 
form, 



P = P(p,e mt ) . 
e tot is given by the equation for the total energy, 



pe tot = pe mt +p- 



222 
v l +v 2 + v 3 



(16) 



(17) 



where Vu U2, v 3 are the components of the velocity vector, and <t> is the gravitational potential. In C05BOLD, a 
prescribed, time-independent gravitational potential is used, so far. Self-gravity is not accounted for. The gravity field 
is given by 







( j_ ) 












d 


92 






, 93 ) 




d 







(18) 



With C05BOLD, Eqs. (TT~3T> - cTT~5T> are solved with the hydrodynamics module described in Sect. 13.51 
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The equations of ideal magnetohydrodynamics (MHD), including gravity and radiative energy exchange, are 
written in the more compact vector notation as 

dp 

■£ + v ■ (pv) = o , 

at 

+ V-(pvv + (p+ B — ^I-Bb) =pg, 



dt V V 2 

<9B 



dt 



(19) 

+ V ■ (vB - Bv) = , 



3gS + V ■ ((p etot + P + 5-5) v - (v • B)B + F rad ) = . 

Here, B is the magnetic field vector, where we have chosen the units such that the magnetic permeability p is equal 
to one. I is the identity matrix and a ■ b = 2a a kbk the scalar product of the two vectors a and b. The dyadic tensor 
product of two vectors a and b is the tensor ab = C with elements c mn = a m b„ and the nth component of the divergence 
of the tensor C is (V • C)„ = Ytm dc m „/dx m . In this case, the total energy is given by 

v v B B 

pe tot = pe m +P—;- + — — + P® > (2°) 

where e mt is again the internal energy per unit mass. The additional solenoidality constraint, 

V-B = , (21) 

must also be fulfilled. The equation of state and the equation for the gravitational field are given by Eq. ( TToT l and 
Eq. ( [T8T l. respectively. With C05BOLD, the equation system, Eq. ( fT9l ), is solved with the MHD module described in 
Sect.lTTl 

In addition, there are equations for the non-local radiation transport solved with C05BOLD with the modules 
described in Sect. 13.6.31 and Sect. 13.6.41 These modules account for the frequency dependence of the opacities by the 
multi-group technique described in Sect. l3.6.2l In the following equations, the subscript v refers to the index of the 
frequency group. 

The variation of the intensity I v along a ray with direction n can be described by the radiative transfer equation 

— n- Vly + Iy = Sy . (22) 

PK V 

The group-averaged opacities k v are typically given as a function of temperature T and gas pressure P, 

Ky=Ky(T,P) , (23) 

and the group-integrated source function, S V (T), is normalized such that 

V S v = B(T) = -T 4 , (24) 

V 

where B(T) is the frequency-integrated Planck function. Introducing the optical depth t v according to 

dr v = pK v n ■ dx , (25) 
where n ■ dx is the path increment along the ray, the radiative transfer equation can be written as 

^+I V = S V . (26) 
dr v 

The frequency-integrated radiative energy flux vector in direction ft is given by angular integration over the full sphere, 
and summation over frequency groups 
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F r ad = 2j M«)ndfl • (27) 



The energy change due to radiative transfer can then be computed from the flux divergence as 



2rad - -V ■ F ra d 



(28) 



To include additional physics such as chemical reactions (Sect. l3.8Tb , dynamic hydrogen ionization (Sect. l3.8T2t 
or dust (Sect. 13. 831 ) the above equations are augmented by 



where the number densities n, represent the densities of chemical species, ionization states, or dust particles. The 
source term _> , accounts for chemical reactions, ionization and recombination, or dust formation. 

2.3. Basic numerics 

The numerical simulations described here are performed with C05BOLD (Conservative COde for the Compu- 
tation of COmpressible COnvection in a BOx of L Dimensions, L=2,3). It uses operator splitting 19111 to separate the 
various (usually explicit) operators: the hydrodynamics CSect. l3.5l > or magnetohydrodynamics (Sect. 13.71 ). the tensor 
viscosity f Sect. l3.5T6t . the radiation transport (different for the two setups, see below; local models: Sect. 13.6.31 or 
global models: Sect. 13.6.41 ). and optional source steps (e.g., due to time-dependent dust formation or hydrogen ion- 
ization, Sect. 13.81 ). The tabulated equation of state accounts for the partial ionization of hydrogen and helium and a 
representative metal (Sect. l3.4l ). The opacities can be either gray or can account for the frequency dependence via an 
opacity-binning scheme (Sect. 13.631 ). Parallelization is done with OpenMP. 

C05BOLD is used for two different types of model geometries, which are characterized by different gravitational 
potentials, boundary conditions, and modules for the radiation transport: in the local-box (or box-in-a-star) setup 
(Sect. l3.2TI ). used to model small patches of a stellar surface, the gravitation is constant, the lateral boundaries are 
periodic, and the radiation transport module relies on a Feautrier scheme applied to a system of long rays (Sect. 13. 631 ). 
In contrast, supergiant simulations employ the global or star-in-a-box setup (Sect. l3.2T2l for which the computational 
domain is a cube, and the grid is equidistant in all directions. All outer boundaries are open for matter and radiation. 
The prescribed gravitational potential is spherical. For this setup, a different radiation-transport module is used, which 
implements a short-characteristics method (Sect. 13. 631 ). 

Some more technical informations can be found in the C05BOLD Online User Manual Fl 

3. Detailed numerics 

In this section, we present some numerical details of the code that are adapted to the conditions found in stellar 
atmospheres. 

3.1. Numerical grid and independent variables 

Instead of the conserved quantities, Eq. (TT2l . we choose the primitive variables 



as independent quantities, using integer indices for the components of a vector. Since the conserved variables are 
purely algebraic combinations of the primitive variables, the primitive variables can be directly updated using the 
conservation laws Eqs. ([T3l)-(fT5b or Eqs. ([T9T > without dismissing conservation-law principles. This is explained in 
more details in Sects. [3331 and l3~5. 41 

The hydrodynamics variables p,v\,V2,v 3 , and e; nt , are cell centered with grid coordinates (;c i,x C 2, x C 3), whereas 
B\,B2, and B3 are cell-boundary centered with coordinates (x\,\,x\q, x\a). The grid is Cartesian. The grid spacing 
may be non-equidistant. Additional subscripts are used to describe the grid indices. The hydrodynamics variables 
p, v\, V2, U3, and e m must be thought of as cell-averaged quantities, while #2, and B3 are mean magnetic flux den- 
sities through cell interfaces. 



See http://www.astro.uu.se/~btyco5bold_main.html http://www.co5bold.com 



^ + V.(„,v) = S ; . , 



(29) 



p, v u v 2 ,V3, eintG B U B 2 ,B 3 ) 



(30) 
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3.2. Boundary conditions and setup 

Global models, that simulate an entire star-in-a-box (typically a red supergiant, Sect. l4.7t . differ essentially in 
boundary conditions and the gravitational potential from local box-in-a-star models, that simulate only a small piece 
of a star close to the main sequence. The fundamental parameters are the effective temperature, T e ff, describing the 
radiative flux per area in local models, or the luminosity in global models, the surface gravity, g, and the chemical 
composition of the stellar material. 

3.2.1. Local models 

Local box-in-a-star models are designed to simulate a small patch at the surface of a star, ignoring effects of the 
spherical geometry and variations in gravity. The computational domain is a Cartesian box with constant, downwardly 
directed gravitational acceleration given by 

g = (0,0,-0) ■ (31) 

The side boundaries are usually periodic. Closed walls are a rarely used option, as they tend to attract downdrafts. 
The top boundary is generally either hit under some finite angle by an outgoing shock wave or it lets material 
fall back into the computational domain (often with supersonic velocities): there is not much point in tuning the 



formulation for an optimum transmission of small-amplitude waves 19211 . Instead, a simple and stable prescription 
that lets the shocks pass is sufficient. It is implemented by assigning typically two or more layers of ghost cells 
(the number depending on the order of the reconstruction scheme), with boundary values, for which the velocity 
components and the internal energy are kept constant. The density is assumed to decrease exponentially with height 
in the ghost layers, with a scale height set to a controllable fraction of the local hydrostatic pressure scale height. 
The layers of ghost cells are located outside the computational domain proper. The control parameter allows for the 
adjustment of the mean mass flux through the open top boundary. 

The bottom boundary of a standard solar model is located well inside the convection zone, where the material 
coming from below is assumed to have the entropy of the adiabat of the deeper convective envelope l3~5ll . The 
corresponding boundary condition prescribes the entropy of the ascending material, ensures a zero total mass flux, 
and reduces pressure fluctuations for stability reasons. Horizontal velocities are assumed to be constant with depth. 
The values of p, e m , and the vertical velocity U3 in the lowermost grid layer are actually modified during the application 
of this boundary condition. Therefore, the conservation laws are only valid in the volume above the bottom layer. For 
each cell in the bottom layer the following steps are performed: 
The equation of state is solved, 

EOS(p,e int ) -» s, P, T, r u r 3 , c s , (32) 

to get the entropy, pressure, temperature, first and third adiabatic coefficient, and the sound speed. Horizontal averages 
of the density and pressure (p) (0) , (P) over the entire bottom layer are computed, where the superscripts (0) , . . . , <3) 
here and in the following equations denote the sub step. A characteristic time scale is estimated by 

(t cbdX ) = Ax 3 /(c s + N> ■ (33) 

In cells with an upflow (1)3 > 0), mass and energy are modified according to 

(1) At -p 2 t (r 3 - l) 

p y = p + Change" ^ (inflow ~ S) , (34) 

'char " 1 1 

efl = Cint + Change" T (l ^ I (^inflow ~ S) (35) 



?char \ r\ / 

with the two external parameters Change (~0.1) and inflow The latter controls the effective temperature T e ff. To 
reduce deviations of the pressure from the horizontal mean, the following corrections are applied to all cells in the 
bottom layer: 

p <2> = p (i) + Crehange I «P> - P) , (36) 

'char C s 
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4 2 t = 4t + Cpc^ — j—dP) - P) , (37) 

adding another parameter exchange (~0.3). To keep the total mass in the model volume unaltered, the density in the 
bottom layer is corrected with 

P 0) =P (2) +(p) m -(p {2) ) ■ (38) 

Because of this step, this boundary condition acts as a closed boundary for plane-parallel waves. Finally, the vertical 
velocity is modified to ensure a zero-average vertical mass flux, 

(1) (P (3) /QQ\ 

Vi = V3 - IpW ■ m 

Now, the old values are replaced by the new ones, 

(new) = (3) (new) = (2) (new, = (1) (4Q) 
" " ' int int ' 3 3 v / 

Later, during the hydrodynamics step, the ghost cells are simply filled with constantly extrapolated values from the 
bottom layer while keeping the gravitational potential constant in these layers. 

3.2.2. Global models 

For global models, the gravitational potential depends on the radius r only. The 1/r potential is a good approxima- 
tion for the outer layers of supergiant stars, which have a small massive core surrounded by an extended low-density 
envelope. To avoid the central singularity the potential is smoothed near the center. The potential can also be flattened 
at large distances to artificially enlarge the pressure (and density) scale height preventing extremely low pressures and 
densities in the corners of the simulation box. The potential is given by 

0>(r) = -GM„(^ +r 4 / Vl +('-M) 8 ) _1/4 , (41) 

where M» is the mass of the star to be modeled and ro and r\ are smoothing parameters in the core and the outer 
envelope, respectively. Within the sphere r < Tq, a source term to the internal energy provides the stellar luminosity. 
Motions in the core are damped by a drag force to suppress dipolar oscillations. 

All six surfaces of the computational box employ the same open boundary condition, which is also used for the 
top boundary in the local models ("Sect. l3.2TTb . 

For global models the temperature/pressure range of the photospheric opacity tables is insufficient. It is therefore 
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merged at around 12 000 K from high-temperature OPAL data 19311 and low-temperature PHOENIX data 19411 



3.3. Initial conditions 



Due to the chaotic nature of stellar convection 09511 and the primary interest in averaged or statistical properties, the 
details of the initial conditions hardly matter, except for initial strong magnetic field configurations. On the other hand, 
the total mass within the computational domain is of main importance. However, choosing a pressure and temperature 
distribution too far off from the (usually close to hydrostatic) mean conditions requires an unnecessarily long time, 
until plane-parallel pulsations have settled down and the stratification is thermally relaxed. It is often advisable to 
start with a standard ID atmosphere model (e.g., produced with PHOENIX as in I960), to expand it trivially into the 
second and third dimension and to add small velocity fluctuations to it as seed for convective motions. An even better 
alternative is to use an existing 3D snapshot with similar parameters - if available - and scale it to the desired model 
properties. 

Even with a careful construction of the start model, transient plane-parallel pulsations are common. These pulsa- 
tions are generated by tiny deviations from the exact numerical hydrostatic equilibrium in the deeper layers, causing 
larger amplitudes in the tenuous top layers. To damp them out, a vertical drag force acting only on the horizontal 
average of the vertical mass flux can be applied in the initial phase of a simulation. 



13 



3.4. Equation of state 

Under the conditions of cool stellar surfaces, a lot of energy can go into the ionization of hydrogen and helium. 
In C05BOLD, the equation of state (EOS) accounts for the ionization balance of HI, HII, H2, Hel, Hell, Helll, and 
a representative metal. Pre-tabulated values as functions of density and internal energy are used (log p, log e mt — > 
logP, log T, s). In fact, the coefficients for a bicubic interpolation of (log P, log T, s) are stored. Thermodynamic 
derivatives are computed from the corresponding derivatives of the polynomials. 

3.5. Hydrodynamics 

In general, a hydrodynamics scheme should 

1 . be consistent with the original hydrodynamics equations, 

2. be stable, 

3. solve the hydrodynamics equations in 3D with reasonable accuracy, i.e., be of high order whenever possible 
and represent discontinuities with only a few grid points, 

4. be conservative to handle shocks properly and give constant total fluxes in stationary cases, which is particularly 
important for modeling convection, 

5. include source terms due to gravity in a proper way to allow static solutions, so that especially the construction 
of an exactly hydrostatic stratification in radiative equilibrium is possible, 

6. handle a general equation of state (from a table), 

7. be fast, e.g., easy to vectorize, to parallelize, and to make proper use of the various CPU caches, 

8. handle various geometries (in this case ID, 2D, and 3D models), 

9. be not too complex but stay fairly simple, 

10. allow the coupling with additional physics (especially radiation transport). 

Solvers differ in how close they get to the individual design goals. For instance, total energy conservation might 
get sacrificed to improve the code stability in cases of large Mach numbers. And with detailed (read, time consum- 
ing) radiation transport modules, the performance of the (usually comparably fast) hydrodynamics modules becomes 
unimportant. 

The hydrodynamics scheme of C05BOLD uses a finite-volume approach. By means of operator (directional) 
splitting 09 ill , the 2D or 3D problem is reduced to one dimension. To compute the fluxes across each cell boundary in 



every ID column in x c \ direction, an approximate ID Riemann solver of Roe type 19711 is applied, modified to account 
for a realistic equation of state (Sect. l3.53T l. a non-equidistant grid (Sect. l3.5TTT i. and the presence of source terms 
due to an external gravity field (Sect. l3.53l l. The partial waves are reconstructed and advected with upwind-centered 
fluxes. A slope limiter (MinMod, SuperBee, but usually van Leer) ll98ll or a reconstruction with monotonic parabolae 
(Colella and Woodward Isstl) is applied to decrease the order of the scheme in the neighborhood of discontinuities for 
keeping it stable while preserving higher-order accuracy in the case of smooth flows. 

The standard Roe solver has been extended in several ways to fit the particular problem of stellar surface convec- 
tion as is explained in the following subsections. 



3.5.1. Non-equidistant grid 

The hydrodynamics scheme handles Cartesian grids only. They may be non-equidistant in any direction. Without 
gravitation, the location of the cell centers x c \ has not much relevance as all quantities are either integral values within 
a cell (for instance the mass density) or located at the cell boundaries x\,\ (for instance the mass flux). In this simple 
case, a non-equidistant grid would only have an effect on the reconstruction equations. 

With the inclusion of gravity however, the potential energy within each cell is located at x c \. This means, that the 
pressure should also be located there in order to allow for a correct balance of acceleration due to the pressure gradient 
and the gradient of the gravitational potential. 

The relative position of jc c i t ; within Xbi,; and Xb\j+\ is not set by the hydrodynamics scheme and can be chosen 
within reasonable limits according to the requirements, e.g., of the radiation-transport scheme. 
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3.5.2. Update of the mass density and the velocity 

Given the update of the density in one coordinate direction for cell i in conservation-law form, 

Ap, = -£ (/ P ,,+i - /p,,) , pf ew) = pf d) + Ap, , (42) 

where f P) , is the mass flux in this direction between cell i— 1 and cell i, the update for the momentum in the 1 -direction 
can be reformulated in terms of the update for the velocity as follows. From the conservation-law form 

(pu^ = (pv^ - £ (f pVuM - f pvui ) + AtS pVui , (43) 

where f pVl j is the 1-momentum flux in the considered direction and S pvi j the source term for the 1-momentum, we 
obtain 



(new) (old) 
Vi) '=Vl) 



^ (/p Ul ,i+l - fpv u i) ~ Af S puui + Api Vif U) 



(new) 

Pi 



(44) 



Ap, and p* new) on the right hand side of Eq. d44"l) are known from Eq. d42"l) . Then, the momentum {pv\ ),- = p; v\ , is, up to 
the source term, a strictly conserved quantity. The fluxes f p and f pUl are defined at the cell interfaces and determined 
by an approximate solution of the Riemann problem as explained in Sect. I3.5l and Sect. 13.5.51 The advantage of his 
formulation becomes apparent when treating the gravitational potential in the derivation of the discrete equation for 
the internal energy. This is explained in Sect. 13.5.41 

3.5.3. Gravity 

The gravitational source term in Eq. (TT4b and Eq. fl9[ destroys the hyperbolic character of the corresponding 
system of equations and inhibits the direct application of an (approximate) Riemann solver. On the other hand, the 
separation via operator splitting is not a good idea in this case, because in stratified atmospheres the pressure gradient 
and the gravity tend to cancel each other (nearly). Their application in sequence - and not together in a single step 
- would cause spurious unwanted accelerations back and forth. On the other hand, the naive combination of the Roe 
solver with the source terms due to gravity into a single operator by simple addition, leads to problems because the 
Roe solver interpretes the strong pressure gradient in a stratified atmosphere as indication of a shock wave, which is 
then treated as such, causing spurious - possibly large - velocity fields. 

There is some freedom in the choice of the exact reconstruction of quantities inside the cells (Mellema et al. llOdO . 
which is used to amalgamate the hydrodynamics with the gravity operator by reducing the pressure jump across a cell 
boundary to the deviation from hydrostatic stratification. The latter is subtracted from the actual pressure inside the 
cells during the computation of the amplitudes of the partial waves. The idea is, that only pressure deviations from 
hydrostatic equilibrium - and not just pressure gradients - should give rise to fluxes of the partial waves. In the exact 
hydrostatic case, the Roe solver should "see" no sound waves. This construction does not supersede the usual source 
terms due to gravitation. 

3.5.4. Update of the internal energy 

The discrete form of Eq. (TT3T > for the total energy is given in conservation law form by 



peint +P z +P®c 



x (new) 



peint +P Z +P$c 



(old) 



- £ + /o,, + i) - (fa + /*,;)) . (45) 



where f,j is the ID flux of the total energy without the potential energy from cell i— 1 into cell i provided by the Roe 
solver, /o,, is the flux of the potential energy and Ci , is the gravity potential in the center of the cell. Using Eq. d42l 
and defining the flux of the potential energy as 

Mi = ®b,if P ,i , (46) 



15 



where <t> D; , is the gravity potential at the interface between cell i— 1 and cell i, all terms containing the gravity potential 
can be combined to yield 



(pem)) 



(new) 



/ \(°ld) , 



-— [(<h.i+l - ®c,i)fp,M ~ (<&b,i ~ ®c,i)fpj] ■ 
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(47) 



The presence of the new kinetic energy at the right side of Eq. d47b does not make the scheme implicit, since the 
velocity does not depend on e ml - it is known from Eq. (144-b and corresponding equations for i>2 and U3. We note that 
this update still conserves the total energy, Eq. ( fTTb . to machine accuracy. The conservative inclusion of the radiative 
energy flux into the energy equation is treated in Sect. 13.6.31 



3.5.5. General equation of state including ionization 

Several extensions of the Roe scheme for a general equation of state have been proposed, see 1 101] and references 
therein. The differences compared to the case of an ideal gas with constant y manifest themselves in the need of 
additional Roe averages depending on which variables are used in the equation of state. For a general equation of 
state of the form Eq. ( TToT ). averages of p, e; nt and of the pressure derivatives and ^-j are required to build the 

Roe matrix. Choosing the usual Roe average for <?j nt and the average yjpip r for the density, where pi and p r is the 
density on either side of the cell interface, the condition 



AP 



dP 

deint 



p op 



Ap 



(48) 



must be fulfilled by the averages of 



and 



where AP, Ap and Ae mt are the jumps of the pressure, density, 



and internal energy at the cell interface. Eq. d48b is not sufficient to determine the averages of the pressure derivatives. 



JUL 

de-m 



Glaister lll02ll suggested formulas for |£ | and 1 which meet Eq. (l48t exactly. However, Glaister's formulae 



<p 



do not lead to the same averaged sound speed as the original formulae of Roe in the simple case of constant y. 
Furthermore, they may also produce unphysical average states llOlll . In C05BOLD, the averaging of the pressure 
derivatives is avoided. Instead, the dimensionless quantities Fi and 3X0 averaged with the usual Roe weights, 

which ensures the consistency with the simple gas case. 

The pressure derivatives in the Roe matrix lead to an additional term 



/ dP 


dP 




\ e ^- p Tp 




;') 



|a (5) 



(49) 



in the energy flux. This term, which vanishes in the perfect gas case, is treated as a contribution from a sixth partial 
wave a (6 \ The wave strength of the entropy wave or (5) is given by 



Ap- 



AP 



(50) 



The term in Eq. ( |49l can become small, possibly causing numerical errors. Using thermodynamic relations, a (6) 
can be transformed into 



a (6) = A(pe mt ) - 



dpe n 



BP 



AP . 



(51) 



which uses a better behaved derivative. 
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3.5.6. Tensor viscosity 

In addition to the stabilizing mechanism inherent in an upwind scheme with monotonic reconstruction, a 2D or 
3D tensor viscosity can be activated. It eli minat es certain errors of Godunov-type methods occurring in the case of 
strong velocity fields aligned with the grid II 1 0311 . Other types of problems can occur when e.g., a shock, which has 
a strength that could easily be handled by the hydrodynamics scheme alone, gives rise to so large opacity variations 
that the radiation-transport routines might get unstable (Sect. l3.6l l. 

To overcome such (possible) problems, an additional tensor-viscosity sub step was included in the code, that can 
add dissipation in a way the Roe solver by its own is not able to produce. The kinematic viscosity is 

v = - (Axi + Ax2 + AX3) Ciinear c s + min (Axi, Ajc2, AX3) max(Axi, Ax2, AX3) (52) 
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with the parameters Cii ne ai-, Canificiai, an d Cs maeon nskv for the linear, artificial (von Neumann type) viscosity, and turbu- 
lent subgrid-scale viscosity (Smagorinsky 110411 '). respectively. Typical values are (0, 0.5, 0.5), respectively. Models 
of solar granulation and similar easy cases do not require this extra viscosity. However, it is usually activated to avoid 
the necessity to tune the numerical parameters individually for each stellar model. For instance the models of the 
more dynamic atmospheres of red supergiants (Sect. l4.7l i need some amount of extra dissipation provided by the ten- 
sor viscosity. Note that the tensor viscosity should not be mistaken for a hyperviscosity. The task of hyperviscosities 
is in C05BOLD done by the reconstruction schemes (MinMod, SuperBee, van Leer, PR etc.). 



3.6. Radiation transport 
3.6.1. Introduction 

In dynamical simulations which take time dependence and coupling between radiative energy transfer and hydro- 
dynamics equations into account, the emitted intensity is only a by-product. Important is instead the energy change 
per numerical grid cell due to the difference of radiative gains and losses. The requirement to solve the radiation- 
transport equations for many grid points and many time steps calls for severe simplifications as, e.g., the restriction 
to gray opacities or to a few frequency groups CSect. l3.6T2b or the treatment of scattering as true absorption. Actually, 
this can make code development easier. But still, there are additional demands on the algorithm: the scheme should 
conserve the total energy, i.e., internal sources and sinks of energy minus losses through the surface should exactly 
sum up to zero. The scheme has to be stable enough to handle complex structures which may sometimes be poorly 
resolved (e.g., chromospheric shocks, Sect. l4.5b . 

Some cases pose only low demands on the complexity of the algorithm: if the entire model is optically thick, a 
diffusion approximation using only differences between neighbor cells is adequate to compute the radiative flux. This 
results in a stable scheme, if the time step is properly limited. If the whole numerical domain is optically thin and the 
radiation field is simple enough, a local cooling function might be sufficient to model radiative energy losses, calling 
for a scheme that is stable if the time step is small enough. 

However, stellar atmospheres are per definition at the transition between optically thick and thin regions. The 
main form of energy transport switches from convective plus radiative in the interior to mainly radiative in the outer 
layers, where mechanical energy fluxes become very small due to the low material density. Still, mechanical energy 
fluxes might be sufficiently large to affect the temperature structure of the chromosphere (Sect. [43), for example. 
While radiative energy transport in the stellar interior can be properly described by local physical quantities through 
the diffusion approximation, in the outer layers, radiative energy exchange occurs non-locally. This means that the 
local radiative flux depends on the physical state of the material in the wider surroundings of a size depending on the 
mean free path of the photons. Large opacity variations due to changes in the ionization states of major constituents 
or due to shock waves can cause changes in the source function on small spatial scales. In numerical models, these 
two effects, amplified by fluctuations in heat capacity, can cause enormous jumps in the radiative relaxation time scale 
from grid point to grid point. 
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Even if a standard scheme is able to overcome all these difficulties and to compute accurate intensities for given 
opacities and source function, there are still several possibilities how to derive the induced energy change per cell, as 
detailed in the following: 

The energy flux through the cell boundaries can be computed from the intensity field, which then gives the di- 
vergence of the flux for each cell, and hence the energy change per cell according to Eq. d2"8l . This would guarantee 
the conservativity of the scheme. But unfortunately, the latter step requires an extreme numerical precision in opti- 
cally very thin regions where the relative flux changes from cell to cell are tiny. A high accuracy is also necessary 
in optically very thick regions where only small deviations of the intensity from the local source function contribute 
to the net flux. Calculating a discrete derivative naturally amplifies noise and has centering problems, when e.g., the 
intensity field is given at cell centers and the divergence is required at the same position. 

Another possibility consists in deriving the energy change per cell from the difference of the angular mean of the 
intensity, 

I /„ dn = — I 

4?r J 4n 

and the local source function S v . Using Eqs. d22l i. d27l i. and d28l l. one obtains 



= TZ f 7 " dQ = ^ J J 7 " sin 9 dd d f ' ( 53 ) 



Qmi = Yj 4m yP Cv-Sy) ■ (54) 

V 

This scheme does not have an explicit conservation form and in fact it will most likely not be conservative. This 
happens because the distribution of the source function within the cell which is used in the integration process for the 
intensity (where typically some high-order interpolation of the source function with optical depth is performed) is not 
exactly the same distribution as is used in computing the difference J v - S v (where the source function is assumed 
to be constant). Another problem of this scheme is the accuracy in optically very thick regions, where numerical 
cancellation may occur between J v and S v . A more indirect way is to derive from the intensity field some geometrical 
information about the radiation field in the form of Eddington factors. These are inserted into the eq uatio ns describing 
the radiation transport v ia th e Eddi ngton moments (see, e.g., for two dimensions Stone et al. 110511 and for one 



dimension Hofner et al. 110611 . 1110710 . This method requires a non-trivial solver to get the radiation field in the first 
place and later an algorithm to solve the huge system of Eddington equations. This procedure might not suffer from 
the problems mentioned above. However, it seems somewhat inefficient first to compute the intensity distribution of 
the radiation field in some detail and then to throw away most of the information and retain the Eddington factors 
only, which are used to solve the radiative transfer equation again - just in a different form. The extra effort can be 
justified by gains du e to, e.g., an elegant handling of scattering processes or the achievement of large time steps with 
an implicit operator Il06ll . though. 

In Sect. l3.6.3l and Sect. l3.6.4l we present two radiation-transport schemes implemented in C05BOLD, that over- 
come the aforementioned problems in different ways. They compute the contribution to the energy change per cell 
on-the-fly during the integration of the intensity for each direction. For standard local-box models with periodic 
boundaries (Sect. l3.2TT1 ). we use a long-characteristics scheme, described in Sect. l3.6.3l while for star-in-a-box mod- 
els with all-open boundaries (Sect. l3.2T2l ). we use a short-characteristics method, outlined in Sect. l3.6.4l 

3.6.2. Opacity binning 

The rate of the radiative energy exchange is highly variable over the relevant spectral range, since the absorption 
coefficient strongly varies with frequency due to the presence of spectral lines, on top of the more gradual change of 
the continuous opacity. In cool stars like the Sun, spectral lines count in the millions so that an exact treatment of 
the frequency dependence in a complex multi-dimensional geometry is beyond present computer capacities, and one 
has to resort to an approximate treatment. An important simplification stems from the fact that one is not interested 
in the detailed frequency dependence of the heat exchange between stellar plasma and radiation field but only in its 
frequency-integrated total effect, Q m d- Nowadays, all multi-dimensional hydrodynamical stellar atmosphere codes 



employ the so-call ed op acity-binning t echn ique. The method was first laid out by Nordlund [23], and later refined in 



works by Ludwig lll08ll . Ludwig et al. jl09ll . and Vogler jl loll . At present, opacity sampling - a statistical technique 
widely applied in standard ID model atmospheres - is discussed as possible replacement of the opacity binning due 
to its better controlled accuracy and greater flexibility. 
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The basic idea of opacity binning is the classification of frequency points by the similarity of their associated 
A-operator - the operator relating source function S v and mean intensity J v of the radiation field: 



We added an index v to the A-operator to emphasize that its form, written in geometrical coordinates, is different 
for different frequencies due to opacity variations. However, in cases where the operator happens to be similar, its 
linearity allows to operate on the sum of the source functions to obtain the integrated mean intensity, symbolically 
expressed as 



where A is some suitable mean of Ai and A2. The problem now is to classify all frequencies into distinct sets 
f2, grouping together as similar as possible A-operators. The A-operator can be calculated from the monochromatic 
optical depth scale r v so that the classification can be equivalently done by grouping frequencies with a similar relation 
between geometrical and optical depth scales. This is in fact the way how one proceeds in practice. 

When trying to classify the frequencies, one is confronted with the problem that the optical depth scales depend 
on the atmospheric model under consideration, i.e., its geometry, the ensuing thermal conditions, and velocities. One 
has to choose a reference model for which the classification is performed. Naturally, this reference model is chosen 
to be close to the stellar atmosphere to be simulated, in the simplest case a ID model of the atmosphere in question. 
Other choices are possible, but in any case, the resulting classification is optimized for a particular set of atmospheric 
parameters and has to be repeated when numerical simulations in other parameter regimes are conducted. Since even 
for fixed atmospheric parameters a large variety of different thermodynamic conditions are met along various lines- 
of-sight in a numerical model (with correspondingly different t v ), limits to the achievable accuracy by the opacity 
binning have to be expected. Thus only a reasonable similarity among r v -scales within an opacity bin is aimed at in 
practice. Typically one is content if the r v -scales of a group of frequencies share the property to reach unity within a 
given range of depth - usually defined via the frequency-independent Rosseland optical depth. This emphasizes the 
emergent radiation intensity as the primary quantity to be captured correctly, obviously an important quantity linked 
to the overall flux properties of a stellar atmosphere. Each opacity bin defines a corresponding frequency group. 

At present, typically between four and twelve frequency groups £2, are used, depending on the desired precision. 
An estimate of the precision is obtained by comparing the integral radiative heating (or cooling) rates obtained from 
the binned opacities with the result obtained at high frequency resolution, both as a function of depth in the reference 
structure used for defining the opacity bins. The estimate relies on the assumption that the reference structure is 
indeed representative of the conditions encountered in the flow simulation. Some refinements to this basic scheme 
are nowadays often added. For instance, it is sometimes advantageous to split an opacity bin as defined before 
into frequency sub-groups, with the idea to separate frequency points which systematically heat or cool particular 
atmospheric layers. This helps to improve the overall energy exchange budget. 

An example is given in Fig. [4] illustrating the results obtained for the ID solar reference atmosphere. The basic 
5-bin/5-group scheme is clearly superior to the gray approximation. The more sophisticated 9-bin/12-group scheme, 
in which three opacity bins are split into two frequency sub-groups, performs very satisfactory and almost perfectly 
reproduces the "exact" heating rate. 

The binned opacities are obtained from a suitable average of the opacities in a particular frequency group and 
stored in look-up tables as a function of thermodynamic variables - in C05BOLD as a function of gas pressure and 
temperature. In addition, the Planck function (as source function), integrated over the frequencies of a group, is stored 
as a function of temperature. This approach only works if the opacities and the source function can be calculated from 
the thermodynamic conditions alone, i.e., are thermodynamic equilibrium quantities. While this is often fulfilled to 
good approximation, there are exceptions. For instance, the formation of dust clouds in cool stellar atmospheres is a 
non-equilibrium process (Sect. l3.83l l. and actual particle properties are only known after solving the governing kinetic 
equations, taking into account the history of the evolution of a particular mass element in the flow. In C05BOLD, 
we proceed by separating the equilibrium part (gas opacities) from the non-equilibrium part (dust opacities). The gas 
opacities are binned into frequency groups in the usual way, and the dust opacities are calculated during the simulation 
on-the-fly and added to the gas opacities. Obviously, this increases the computational demands. 

All in all, opacity binning has been and still is working perhaps better than one might expect from the numerous 
approximations behind the construction of the scheme. Opacity binning has proved to be an efficient way to include 



/v = Av [S v ] 



(55) 



J 1+2 = Ji +J 2 = M [Si] + A 2 [5 2 ] ~ A [Si +S 2 ] 



(56) 
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Figure 4: Performance of the opacity-binning scheme, illustrated for a ID solar model atmosphere. The net radiative heating rate per unit mass, 
Qmd/p (top), and the bolometric radiative flux, Frad/^o (bottom), and are shown as a function of Rosseland optical depth. In each panel, the results 
from gray (dashed), 5-bin/5-group (dotted), and 9-bin/12-group (solid) radiative transfer are compared with the "exact" solution (diamonds), 
obtained with very high frequency resolution. 



the frequency dependence of the radiative transfer in multi-dimensional simulations. However, as alluded to already 
before, the increased computing power might allow to re-consider the approach trading greater computational costs 
for higher physical fidelity. The path to largest gains needs to be identified yet. 

3.6.3. Long-characteristics radiation transport 

The purpose of this algorithm is to compute the net radiative heating rate per unit volume, Qradfe UjiZk), at the 
center of each cell of the hydrodynamical grid (HD grid). The basic idea is to solve the equation of radiative transfer 
on a system of straight long rays (long-characteristics, LC) running from the upper to the lower model boundary at 
a number of different azimuthal angles <f> and inclinations with respect to the vertical (0 < < n/2). As a result, 
we obtain for each frequency group v and for all bundles of rays with orientation (9, <p) the quantity Q a d y (9, <p) = 
pK v (u v (9, <p) - S Y ) at the mesh points along the rays, where the mean-intensity-like variable u v (9, <p) is the average of 
incoming (/") and outgoing (1+) intensity, u v = (7+ + 1~ )/2 (see Fig.[5]i, S v is the group source function, and p/c v is the 
group opacity averaged over the neighboring mesh points along the ray (see Eq.l63ll. Qradfe Hj, Zk) is then constructed 
by interpolating Q m d, v {9, <p) from the ray system to the cell centers of the hydrodynamics grid, and appropriate angular 
averaging and summation over frequency groups. 

Note that the technique described here basically evaluates Qmd according to Eq. (T54b . It overcomes the difficulties 
explained in the context of Eq. ( f54b by solving the transport equation for the difference between mean intensity and 
source function, p v = u v - S v (see Eq.l57b. which gives accurate values of (u v - S v ) for arbitrarily large optical depth. 
At the same time, it allows Q m d,v to be computed such that energy conservation is enforced (see Eq.l65l). The procedure 
is very similar to that described in l23ll . 
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Figure 5: Schematic illustration of the different grids used with the long-characteristics radiative transfer method. The hydrodynamics equations 
are solved on a Cartesian grid (HD grid, black, dots representing cell centers), while the radiative transfer equation is solved on a system of inclined 
rays (red, dots representing the mesh points used with the Feautrier scheme). The HD grid can be refined in vertical direction by additional j-planes 
(thin, blue) to provide sufficient resolution for strongly inclined rays. The cell centers of the refined HD grid have indices (i, k), the mesh points 
along the rays have indices (m, k). 



To simplify matters, <p is restricted to 0, {\/T)n,n, (3/2)n, i.e., we consider only 2D ray systems in vertical slice s 



along the x and y axis of the hydrodynamical grid. The 6 angles are given by Lobatto's quadrature formula II 1 111 : 
typically, 2—4 non-zero inclination angles are sufficient, in addition to a set of vertical rays. All rays start at the cell 
centers of the uppermost level of the HD grid and follow the specified direction, assuming periodic lateral boundary 
conditions, until they reach the bottom of the computational domain. 

As indicated in Fig. [5] the mesh points along the rays are defined as the intersection points with the z-planes of the 
HD grid. As this recipe would imply a rather coarse sampling along strongly inclined rays, we introduce additional 
horizontal planes such that the geometrical separation of mesh points along the inclined rays remains comparable to 
the vertical resolution of the original HD grid. The coordinates of the ray points, (x^Zk), where m is the ray index 
and k is the depth index of the refined HD grid, are equidistant in x. 

The main steps of the whole procedure may be summarized as follows: first, the source function, S v , and the opac- 
ity per unit volume, p k v , are interpolated from the HD grid to the mesh points of the ray system. Linear interpolation 
of S v and \og{pK v ) is adopted for the vertical direction (additional z-planes), while linear interpolation of S v and p/c v 
is used in horizontal direction. Note that only a ID interpolation along the Cartesian grid lines is required. Give n p k 



on the mesh points along the rays, we represent p k v between two mesh points by a monotonic cubic polynomial [ 112] 
to obtain the optical depth increments At v by analytical integration. Next, we solve the equation of radiative transfer 
along bundles of rays in the form of the second-order differential equation: 

d 2 Pv d 2 S v _ 

~dr[^ Pv ~~dTl ' P v = Uv ~ Sv • 

where t v is measured along the (inclined) rays. This modified Feaut rier e quation is solved by the forward-elimination 
and back-substitution formalism originally described by Feautrier 1 113 1 (see also Mihalas giving p v {x m k,Zk) at 
the mesh points of the ray system. 

At the lower boundary, where conditions are optically very thick in general, we can choose between two basic 
options: if the bottom layer is located in a radiative zone, and we want to enforce a given radiative flux F r ad,v through 
the lower boundary, the condition is 

dp v 3 F mdtV dS v 

— = cos(0) - — . (58) 

dr v 4 n ar v 

If the bottom layer is located in a convective zone, where the radiative flux through the lower boundary is negligible 
compared to the energy flux carried by the flow, a reasonable boundary condition is to require the net radiative energy 
exchange to vanish in each frequency group, 

V • F rdd , v = or p v = . (59) 
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Note that this does not imply F m d,v — 0. 

At the uppermost layer, the optical depth is computed as 

t Vj0 = H T p K v fi , (60) 

where H T is the mean optical depth scale height at the top of the model, {H T Y l = -(dln(pA- v )/dz)o ~ -(dlnr v /dz)o. 
The incident radiation is given by 

/;(r v ,o) = (l-e- T '- )5 y) o + /; , (61) 

where 5 Vj o is the mean source function of the upper layer, and /* denotes the incident intensity due to an arbitrary 
external source (usually zero). In terms of p v , the upper boundary condition for radiation can be formulated as 

Pv - ^ = (1 - e- T *>) 5 V , + / v * - S v + ^ . (62) 
dr v dr v 

Next, the quantity q v is computed at all mesh points of the ray system as 

. . „ T v (x m- k+l, Zk+l) - Tv(Xm,k-l, Zk-l) , , 

q V {x,nk,Zk) = COS0 ■ Pv{x mk ,Zk) 

Zk+l ~ Zjfc-1 

= {pl^{U V ~ Sy)} mk . (63) 

Finally, the partial heating rates q v are interpolated back onto the HD grid in a conservative way, such that for all 
height levels k 

^ qv(xmk,Zk) = ^ q v (xi, Zk) ■ (64) 

m i 

Qmd(xi, yj, Zk) is then built up by adding the individual contributions q v (xi, yj, Zk) of the different ray directions (8, (f>) 
with their appropriate integration weights, and summation over all frequency groups v. 

By virtue of the definition of q according to Eq. d63l . and the requirement of a conservative back interpolation as 
expressed by Eq. (l64l . our long-characteristics radiative-transfer scheme conserves energy in the sense that for each 
frequency group 

f f F l Z v dxdy- [ [ F^ v dxdy = f f f Q 1 - ad , v (x,y,z)dxdydz . (65) 

Here, v and are, respectively, the net radiative energy flux through the upper and lower boundaries of the 
model, computed directly from the ray system intensities at the top and bottom level. Note that Eq. d65l l holds only if 
the volume integral includes the Q ld d,v obtained at the additional horizontal sub-levels introduced for grid refinement. 
The final Q ml $ on the original HD grid must therefore be computed as a suitable average over the neighboring z 
sub-levels to ensure energy conservation. 

A distinct advantage of the long-ray approach is that it allows an efficient solution of the transfer equation for 
beams of parallel rays by means of the Feautrier scheme, which is very fast and elegant, automatically ensures the 
correct asymptotic diffusion limit at large optical depth, and could easily account for scattering along single rays (for 
an early example of this a pproa c h see Cannon ill 14ll ). In principle, the LC method can also be combined with integral- 



operator techn ique s (e.g. Hi 1511 . Ill 1610 . which, however, are numerically less efficient and suffer from interpolation 
issues ( ill 17ll . 1 1 18]). In contrast to what is assumed in Kunasz and Auer ll 19ll . the computing time of our LC scheme 
scales linearly with the number of HD cells and the number of frequency groups, as for the short-characteristics 
scheme. It scales in a non-linear way with the number of 0-angles, since more-inclined rays are longer and have a larger 
number of mesh points. The computing time can be reduced by computing grad from the diffusion approximation in 
the lower, optically very thick layers of the model. Compared to the ray-system solution, the computation of the 
diffusion approximation comes almost for free. 

A disadvantage of the LC method is the necessity of extensive interpolation from the HD grid onto the ray system 
and back. This procedure is prone to problems with "leaking" of heating or cooling to neighboring cells in the presence 
of localized "hot spots", as described in the following Section [3.6.4l (cf. Fig. [6]). To some degree, such problems may 
be abated, at the expense of higher computational cost, by increasing the number of rays per unit length in horizontal 
direction. 
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Figure 6: Initial run of the source function (top curve), induced changed of energy per cell (center curve), and run of the source function after one 
time step (bottom curve). The three quantities are plotted as function of optical depth for a few grid cells whose boundaries are depicted by vertical 
dotted lines. The values at the cell centers are marked by squares. The two sub-intervals in each cell can have different values of the energy change. 
The lower curve shows the source function after one time step for a constant heat capacity per grid cell (thick line) and the case where the heat 
capacity in the left neighbor of the hot cell is smaller by a factor 10. 




Figure 7: The sketch on the left illustrates the naming convention used for the case of a linear dependence of source function S on optical depth 
T within a single interval as opposed to the case of a piecewise linear source function as in the plot on the right, where the interval [ro.Tj] is split 
into two sub-intervals with width Ato/2 and Ari /2, respectively. The source function varies linearly in each sub-interval. However, it is allowed to 
have a jump at the transition. 



3.6.4. Short-characteristics radiation transport 

The LC scheme described in the previous section is part of C05BOLD since the very beginning. It is adapted to 
the conditions of plane-parallel atmospheres in local models: e.g., it heavily makes use of periodic side boundary con- 
ditions. The angular distribution of rays is chosen to optimize the vertical radiative flux. The diffusion approximation 
used in the deeper layers can save some computational time. 




Figure 8: Short-characteristics step to get the intensity in the target cell (bottom) from the values at the previous cell plane (top) in two dimensions: 
left: standard integration with one short ray and interpolation of intensity and source function at the top grid plane; right: separate rays for each 
neighbor cell with a summing up of the intensity in the bottom grid cell, the splitting of each ray gives the intensity change within each cell. 
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While for local models there was no reason to spend time on experimenting with another radiation-transport 
scheme, this changed for global models where the conditions are different: the vertical direction is not preferred 
anymore. Instead, all sides of the computational domain are open for radiation. The numerical resolutions is in general 
worse than for local models and the violent flows and give rise to large local temperature and opacity fluctuations. 
This means that errors caused by the interpolation in LC schemes would become more apparent. 

The short-characteristics (SC) scheme in C05BOLD overcomes these stability problems at the possible expense 
of the accuracy of the vertical radiative energy flux. The basic idea of not following rays through the entire volume 
is the same as in Kunasz and Auer 11 1911 . But a different way of interpolating the intensity and the source function 
makes it better adapted to the use within an RHD code. 

The main emphasis during the development of the SC scheme has been put on stability by preventing local peaks 
of the source function from "leaking" into neighbor cells and causing an unwanted smearing of the cooling or heating 
term (see Fig. [6}. This requires a special reconstruction of the source function within optically thin cells in the ID 
radiation transport operator and a carefully chosen interpolation within the SC scheme. 

Instead of a Feautrier scheme as in Sect. l3.6.3l the analytic solution of the ID version of the radiation transport 
equation (22\ with linear source function (Fig. [7] left) is used as atomic operator, 

h = h e- AT + S i - So e- AT + ^ [e- AT - l] , (66) 

which guarantees the positivity of the source function everywhere. 

The energy change has to be computed accurately in the optically very thick (e.g., in the center of a toy stellar 
model with At > 10 8 ) and in the optically very thin (e.g., in some regions far away from the surface of a red supergiant 
model with At < 1CT 20 ). Both cases pose no problem for the formal solution because in the former case the intensity 
is essentially given by the local source function. And in the latter case, the changes to the radiation field due the 
contribution of the extremely thin regions can be safely ignored - or simply added to the much larger intensity along a 
ray and therefore absorbed by the limited machine precision. However, optically very thick or thin regions still interact 
with the radiation field and the local heating or cooling is significant and has to be computed in a time-dependent code. 
The SC scheme in C05BOLD uses different arrangements of terms in optically thin and thick regimes to account for 
round-off errors, giving accurate values for optically very thick or thin regions - even running only in single precision. 

Separate integration steps are employed for the interval from the cell center to the boundary and from the cell 
boundary to the next cell center (Fig. [7] right) to get the intensity change within each cell, from which the energy 
change per cell is computed. In the optically thin regime, the slope of the source function is reduced (Fig. [71 right) to 
suppress leaking of cooling or heating from one cell to the next (Fig. [6}. Each ray inclination requires an integration of 
the intensity in both directions. For inclined ray directions, there is one pair of intensity-integration steps for each pair 
of neighbor cells to avoid the leaking of cooling or heating associated with the spatial interpolation of sour ce fu nction 
and/or intensity (Fig. [8]). That means, that in contrast to the well-known SC scheme in Kunasz and Auer ill 19ll . even 
for a single direction there might be more than one ray connecting a cell with its neighbors. 

The numerical scheme proceeds as follows: 

At the beginning of each radiation-transport sub step, the temperature T is computed from density p and internal 
energy e mt for every mesh point of an equidistant 3D Cartesian grid. For every frequency group (Sect. l3.6T2t opacity k 
and source function S are calculated by interpolating in precompiled tables. Next, for every ray inclination the optical 
thickness At of each cell is calculated. 

At the beginning of each integration step, the boundary values of the intensity have to be set. For the SC 
scheme, only open boundary conditions (zero infalling intensity) are implemented (for simulations with periodic side 
boundaries, the LC scheme (Sect. 13. 631 is used, instead). 

The integration proceeds then layer by layer along the axis that is closest to the inclined ray direction. For each ray 
direction, the intensity at each cell does not depend on its neighbors within a layer but only on cells in the previously 
computed layer. That means, that the innermost loop in each layer can be efficiently vectorized. The next loop is 
parallelized with OpenMP directives and the outermost loop performs the integration. 

Each complete 3D radiation-transport step includes directions according to the coordinates of the corners of regular 
polyhedrons, which results in equal weights for all rays. After the loop over all inclinations and the loop over all 
frequency groups, the energy change per time is derived from all the accumulated intensity changes and used to 
update the internal energy ej nt in each cell for given time step Af. Here, "conservation of intensity" translates into 
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"conservation of energy" and ensures the conservativity of the radiation-transport update step (except for the losses 
through the outer boundaries). 

From one sub time step to the next, the orientation of the polyhedron can change randomly to give some coverage 
of the entire sphere. However, some simulations are restricted to rays aligned with axes or diagonals resulting in 
a considerable speed-up while loosing some angular resolution. There are several radiative time steps per hydro- 
dynamics time step possible to compensate for the short radiative time scale compared to the hydrodynamic one. 

In cool supergiants close to the Eddington limit, radiation pressure plays an important role in the stellar atmosphere 
and the wind of asymptotic giant branch ( AGB) stars is driven by radiative pressure on dust. With the scheme presented 
above, the three components of the radiative acceleration can easily by computed from the intensity change per cell. 



3.7. MHD 

In C05BOLD, the numerical scheme used for the solution of the equations of magneto-hydrodynamics is quite 
different from the one employed for the case of pure hydrodynamics described in Sect. 13. 51 In the case of solar and 
stellar magnetoconvection, the scheme must be able to deal with highly stratified flows where the plasma-/? (i.e., the 
ratio of the thermal to the magnetic energy density of the plasma) varies over several orders of magnitude. A special 
requirement of MHD calculations is the enforcement of the divergence-free condition V B = for the magnetic field. 
Violating this condition can lead to unphysical forces, which can degrade the solution jl20ll . Several methods have 
been developed to enforce this condition either to roundoff error or a pproximat ely to the order of the scheme. One 



method is to use the eight-wave formulation of the MHD equations 1121L I122H . The additional wave is associated 



with the propagation of magnetic monopoles. In the eight-wave formulation, additional source terms proportional to 
V ■ B appear, i.e., the equations are no longer conservative. Another method uses a cleanup step at the end of each 
time step, removing the errors in V ■ B = 0. This requires the solution of a Poisson equation at each time step. A 
third pos sibil ity, which is used in the MHD module of C05BOLD, is the constrained-transport method of Evans and 
Hawley IU23I1 . It uses a special finite-difference discretization of the induction equation on a staggered grid such that a 
discrete formulation of the divergence-free condition for the magnetic field is maintained to machine accuracy. All of 
these methods can be treat ed as modifications of an underlying base scheme. A detailed comparison of these methods 
can be found in Toth IU24I1 . 



Another difficulty in MHD simulations is to keep the thermal gas pressure positive B125L 112611 . Since the gas 



pressure is a dependent variable when using the conservative form of the MHD equations, it is computed by subtracting 
the potential, the kinetic, and the magnetic energy from the total energy, e tot . When the magnetic energy is much larger 
than the internal energy, i.e., for small values of the plasma-/?, small errors in the total energy can drive the gas pressure 
to negative values. This can be a problem in the solar chromosphere, where values of (3 « 10~ 4 are common, whereas 
the gas pressure dominates in the sub-photospheric layers where /3 is huge. In the MHD module of C05BOLD, 
several provisions are made to avoid a negative gas pressure. To keep the magnetic field solenoidal, C05BOLD uses 
the constrained-transport method in combination with a Godunov-type finite-volume scheme as the base scheme. In 
the following, each component of the scheme is described in more detail. 



3. 7.1. Spatial and temporal discretization 

The spatial discretization of the MHD equations is similar to the hydrodynamic case, i.e., the hydrodynamic 
variables are cell centered. The magnetic fields are located at the cell interfaces. The cell-centered magnetic field 
components, which are required by the Riemann solver of the base scheme, are computed from the magnetic field 
components at the cell interfaces by linear interpolation. Then all cell-centered variables are updated by the base 
scheme. The extension to second order in space is done by linear reconstruction o f the prim itive variables p, v, B, 
P, and pe mi . Second order in time is achiev ed ei ther by a Hancock predictor step 111271 112811 or by a second-order 
TVD-Runge-Kutta time-integration scheme Il29ll . In some situations, where the second-order scheme would result 
in negative gas pressure, the scheme is locally reduced to first order. 



3.7.2. The approximate Riemann solver 

In the hydrodynamic scheme of C05BOLD, a Roe solver is used (Sect. l3.5l l. However, the Roe solver does not 
guarantee positivity of the density and the pressure. This problem, which is also present in the hydrodynamic case 
gets worse for MHD. Whereas in the hydrodynamic case, reducing the time step often helps to overcome the problem, 
in MHD simulations, the problem remains, even if the time step is reduced considerably. 
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It can be shown that the HLL sol ver 1 1 30ll ensures positivity of the gas pressure and the density if the exact solution 



of the Riemann problem is positive 113111 . For MHD, this is the case only if there is no jump in the normal component 
of the magnetic field. In ID, the divergence-free condition enforces the normal component of the magnetic field 
to be constant. For multi-dimensional problems however, using cell-centered magnetic fields, jumps in the normal 
component of the m agne tic field occur even if the divergence-free condition is fulfilled in a discrete sense. It was 



shown by Janhunen 1126H that allowing magnetic monopoles, which arise from these jumps, and taking into account 
their contribution to the Lorentz force, an additiona l sou rce term occurs in the induction equation only. Using a 
special discretization of this source term, Janhunen Il26ll demonstrated numerically that the HLL solver for MHD 
always provides positive gas pressures. 

We use the method of Janhunen for the MHD module of C05BOLD. However, it should be noted that this source 
term is only used for the computation of the fluxes by the HLL-solver of the base scheme. For the update of the 
magnetic field by the constrained-transport method, this source term is not used so that the magnetic field stays 
divergence-free. 

3.7.3. The constrained-transport step 



C05BOLD uses the flux-interpolated constrained-transport method of Balsara and Spicer 1113211 . First, the electric 
field at the cell edges is computed from the fluxes at the centers of the cell interfaces, provided by the base scheme. 
The magnetic field at the cell interfaces is then updated with this electric field, applying Stokes theorem to every face 
of a cell. The updated cell-centered magnetic field from the base scheme is discarded. The new cell-centered magnetic 
field is computed from the updated magnetic field at the cell interfaces by linear interpolation. 

Since the new cell-centered magnetic field is different from the magnetic field provided by the base scheme, the 
internal energy, e mt , must be modified after the constrained-transport step according to 

B* B - B B 

eint = e int + Y p ' (67) 

where B* and e* are the magnetic field and the internal energy provided by the base scheme. If t his c orrection 



would result in negative gas pressure, it is not performed, i.e., e ml — e* mt (see also Balsara and Spicer 113211 ') and the 
total-energy conservation is sacrificed in favor of improved robustness. 

3. 7.4. Dual-energy method and Alfven-speed reduction 

Even if a scheme guarantees positivity of the gas pressure, this does not necessarily mean that the gas pressure 
is computed accurately. In fact, by using the total energy equation for the computation of the internal energy, the 
discretization errors in the total energy, the kinetic energy, and the magnetic energy of the scheme tend to be imposed 
on the internal energy. One could use the entropy equation or the equation for the thermal energy itself, instead. 
Another possibility is to use the total energy equation in combination with one of these equations. For example Balsara 
and Spicer 112511 use the entropy equation for the update of the internal energy in regions with strong magnetic fields. 
For the MHD module of C05BOLD, the so-called dual-energy method, i.e., a combination of the equation for the 
total energy and the equation for the thermal energy is used. In regions with a large /3, the internal energy is updated 
with the equation of the total energy. In turn, when /3 is small (J3 ;$ 10 -3 ), the equation for the internal energy is used 
at the expense of strict energy conservation. Since typically f3 is small in very restricted regions of the computational 
box only, conservation of total energy is still maintained in most parts of the computational domain. 

In order to avoid extremely small time steps due to the CFL condition when the Alfven speed is high, the Alfven 
speed can be limited by artificially reducing the strength of the Lorentz force by a factor 

2 

/ = - °J£f , (68) 
b; + in 

A Amax 

where l>a is the actual Alfven spee d and UAmax is the desired upper limit of the Alfven speed. The method is similar 



to that used by Rempel et al. 113311 . Of course, caution is indicated when using this method. Obviously, it can hardly 
be used for the study of magnetoacoustic wave propagation. However, it may be perfectly admissible in situations, 
where the low-/? regime is merely included as a buffer region to the (upper) boundary of the physical domain. 
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3.7.5. Ohmic diffusivity 

While it is not necessary for stability, the MHD-scheme of C05BOLD can also handle explicit magnetic diffusion. 
It is treated explicitly in the scheme by modifying the electric field in the constrained-transport-step. A constant mag- 
netic diffusivity and the artificial magnetic diffusivity according to Stone and Pringle ill 3411 are currently implemented 
in C05BOLD. The constant magnetic diffusivity can be used to specify the magnetic Reynolds number 



where r\ is the magnetic diffusivity, L is a typical length scale and u is a typical velocity of the flow. The artificial 
magnetic diffusivity is given by 



where Ax is the grid spacing, j = V x B is the current density, and C is a dimensionless parameter. 

3.7.6. Magnetic boundary conditions 

The boundary conditions for the magnetic field can be specified independently from the hydrodynamic settings for 
the top and the bottom boundaries, and for each of the horizontal directions. Typical horizontal boundary conditions 
used for simulating magnetoconvection in a local box are periodic. Another boundary condition, mostly applied to 
the bottom and the top of the box, consists in setting the magnetic field tangential to the boundary to zero, so that the 
magnetic field lines stay normal to the boundary. A generalization of this boundary condition specifies the obliquity 
of the magnetic field at the boundary. There is also a special condition for the open lower boundary, which allows 
upflows to advect horizontal magnetic field into the computational box. Another boundary condition consists in setting 
the electric field to zero at the boundary. This means that the normal component of the magnetic field at the boundary 
does not change. In this case the magnetic field lines are effectively anchored at the boundary. 

Conditions which keep the magnetic field vertical at the top and bottom boundaries are typically used for the sim- 
ulation of intense, vertically directed magnetic flux tubes in the photosphere of the Sun as they occur in magnetically 
active regions such as plages and enhanced network regions. The advection of weak horizontal field across the bottom 
boundary is used for the simulation of magnetically inactive, very quiet regions on the Sun. With this boundary con- 
dition, it is assumed that convective updrafts transport magnetic fields from deep layers of the convection zone to the 
solar surface. Anchored fields may be useful for anchoring an entire sunspot at the bottom boundary or for anchoring 
horizontal fields at the side boundaries for the simulation of horizontally directed penumbral filaments. 

3.8. Optional modules 

The numerical treatment of the source terms S , in Eq. d29l dealing with different types of dust and chemical- 
reaction networks is implemented as a separate step (see below), following the general concept of operator splitting. 
The optional modules are called after the (magneto)hydrodynamics step for each computational time step. These 
modules treat the mass or number densities of the dust particles or chemical species as additional quantities, which 
are included in the in- and output of the simulation data. Only up to one of these extra modules can be used at a time, 
so far. 

During the (magneto)hydrodynamics solver step, the additional densities are advected with the flow field anal- 
ogously to the gas density. Their transport velocity across each cell boundary is computed from the gas mass flux 
divided by the upwind density, in some cases modified to account for the gravitational settling of dust. 

The contribution of the additional components to the total opacity can be added to the standard equilibrium gas 
opacities (Sect. l3.6T2l >. The boundary conditions are made consistent with the hydrodynamics part of the code. 

3.8.1. Chemical- reaction networks 

Apart from advection across the cell boundaries, the number density «, of a chemical species in a grid cell can be 
changed due to chemical reactions: 
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with the index i for the included species. So far, the implementation is restricted to two- and three-body reactions, 
which is a reasonable assumption for comparatively hot stellar atmospheres. The losses (negative sign) and gains 
(positive sign) due to two-body reactions are described by the first and second right-hand terms with the corresponding 
reaction rates k,2,ij and hi,ju respectively. Analogously, the third and fourth term describe the change due to three- 
body reactions. Such an equation is imposed for each included chemical species, resulting in a system of ordinary 
differential equations of first order. In C05BOLD, the chemical reactions are handled locally for each grid cell 
separately by solving the system of differential equations. It starts with the calculation of the reaction rates k, which 
are functions of the local gas temperature and (for catalytic reactions) also of the number density of a representative 
metal. The influence of the radiation field has been ne glecte d so far. The functions are parametrized with prescribed 



coefficients that are provided in the form of a table [see [133, f° r details] 



The chemical-reaction rates, the number densities of the involved species, and thus their derivatives can differ 
by many orders of magnitude, which can cause the system of equations to be st iff. Thus, an implicit scheme is 
used for the numerical solution. We based our solver on the DVODE package Il36ll with an implicit BDF (backward 
differentiation formula) method and an automatic internal time step. The solution finally provides the number densities 
of the involved species after the overall (global) computational time step. For the numerical simulation of carbon 
monoxide, 7_chemical species and a representative metal are considered, which are connected through 27 chemical 



reactions fl 13511 



Carbon monoxide is a non-negligible opacity source in the solar atmosphere, so that the opacity is in principle 
affected by the deviations from equilibriu m of the CO number density. To account for this effect, t he b ack-coupling 



to the radiative transfer was implemented 113711 . It follows the approach by Steffen and Muchmore H138H . which uses 
two frequency groups. The first comprises the gray Rosseland opacity a-r without the wavelength region around the 
CO fundamental vibration-rotation band in the infrared at a wavelengths around ~ 4.6 fim. This wavelength range is 
simulated with the second band, which is constructed from the gray Rosseland opacity kr and an additional opacity 
Ato- The latter is directly connected to the CO number density that is derived from the preceding solution of the 
chemical-reaction network. 

3.8.2. Tune-dependent hydrogen ionization 

A detailed treatment of the time-dependent ionization of hydrogen is important for atmospheric layers, where 
significant deviations from the ionization equilibrium occur, e.g., in the solar chromosphere (see Sect. l4.51 >. Current 
applications for the Sun use a hydrogen model atom with 5 bound energy levels and a level representing ionized 
hydrogen. The number densities of the individual level populations «, enter as additional quantities in C05BOLD. 
The levels are connected by collisional transitions and 10 radiative transitions (5 bound-bound and 5 bound-free). The 
rate Pjj of a transition between a level i and a level /' is given by Pjj - Cjj + Rjj, where C,y and R,j are the rates due 
to collisional and radiative transitions, respectively. First, these rates are calculated from the local gas density and 
gas temperature, the imposed radiation field, and the level populations and electron densities, that are available from 
the previous time step. Apart from advection, the change of the population number densities and thus the ionization 
degree of hydrogen in a grid cell is then described by a set of time-dependent rate equations of the form 

n, n, 

(S ikon = J] njPji - m J] Pa . (72) 

i*i i*i 

where the terms on the right-hand side are the rates into and out of level ;. In C05BOLD, the set of rate equations 
for all considered energy levels is solved with the DVODE package like that for the chemical-reaction networks (see 
Sect. lSXTY 

An important simplification concerns the usage of fixed radiative rates. In principle, the radiation field and the 
radiative transitions for hydrogen (both bound-bound and bound-free) are connected in such a way that a detailed 
solution has to be found by iteration, which makes it computationally expensive. For the implementation in a multi- 
dimensional radiation MHD code, the rates are fixed and calibrated so that they reproduce the full solution as it is 
implemented in time-dependent ID simulations 

EMEU- 

There is no back-coupling to the equation of state and the 



opacities in the current implementation in C05BOLD, a point, which will be worked on in the future. More details 
are given in Leenaarts and Wedemeyer-Bohm 014111 . 
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3.8.3. Dust 

If the atmospheric temperatures are low enough, not only molecules but larger particles - dust - can form. In the 
Earth's atmosphere, possible types of such particles are for instance aerosols of different compositions, rain, snow, or 
hail, mostly made of water, or particles from volcanic ashes. Many cool stars or substellar objects may contain dust in 
their atmospheres, too (see Fig.fTTb. In the hotter objects, the dust will mostly be made of minerals (e.g., forsterite). 
At lower temperatures (e.g., on the Earth) water plays an important role. For the dust chemistry of the warmer objects, 
it is crucial whether oxygen or carbon is more abundant, because these two elements first form carbon monoxide gas 
(CO) and only the remainder participates in dust formation. The fo rmation, i nteraction, and destruction of grains of 
different chemical composition, size, and shape is difficult to model 1 14211 14311 . Compared to the possible complexity 
of the processes occurring in real objects, the various dust modules in C05BOLD are simple. They are permanently 
under development. Two examples will be outlined in the following: 

Stars at the tip of the asymptotic giant branch lose part of their mass in form of a stellar wind, likely driven 
by radiation pres sure on dust. The formation of carbon-rich dust around such star s was investigated by Freytag 
and Hofner Jlil with C05BOLD and with the 1D-RHD code of Hofner et al. iflQTh . The C05BOLD dust model 
includes a time -dependent description of dust grain growth and evaporation using a method developed by Gail and 
Sedlmayr 1 14511 and Gauger et al. 014611 . In this approach, the dust component is described in terms of four moments 
Kj of the grain-size distribution function, weighted with a power j of the grain radius. The moment Kq represents 
the total number density of grains (integral of the size distribution function over all grain sizes), while the ratio 
K3/K0 is proportional to the average volume of the grains. The equations, which determine the evolution of the dust 
components, are solved considering spherical grains consisting of amorphous carbon. The nucleation, growth, and 
evaporation of grains is assumed to proceed by reactions involving C, C2, C2H, and C2H2. The four moments Kj are 
number densities and are advected with the gas as described in Sect. 13. 81 The gas and dust opacities in this case are 
gray. Some results are shown in Sect. l4.7l and Fig. [TBI 

In contrast to the cool giants, the conditions for the formation of (oxygen-rich) dust in M dwarfs and brown dwarfs 
are fulfilled even in standard ID atmosphere models. However, the comparatively heavy dust grains should sink under 
the influence of gravity and vanish from the visible photosphere, leaving no direct trace in emergent spectra, which 
is at variance with observations. The scheme used in Ir96ll to investigate the question why the d ust does not settle or 
how the material comes back up is based on a simplified version of the dust model used in 810711 . adapted to forsterite 
(Mg2SiO/t, 3.3 g/cm 3 ). In this method, there are only two density fields, one is used to specify the mass density of dust 
particles and the other to describe the m onom ers (the dust constituents), instead of four density fields for the dust and 
none for the monomers as in II 1 0711 and 1 14411 . Therefore, the ratio of the s um of dust and monomer densities to the 
gas density is allowed to change, in contrast to the dust description in lll07ll used for the AGB simulations mentioned 
above. Instead of modeling the nucleation and the detailed evolution of the number of grains, a constant ratio of 
the number of seeds (dust nuclei) to the total number of monomers (in grains or free) per cell is assumed. If all the 
material in a grid cell were to be condensed into dust, the grains would have the maximum radius rd, max , which is set 
to a typical value of 1 /vm. This is close to the typical particle sizes found for the optically t hick part of the cloud deck 
in solar-metallicity brown dwarfs according to the DRI FT-P HOENIX models of Witte et al. 1 14711 . 

Condensation and evaporation are modeled as in 110711 . with parameters and saturation vapor curve adapted to 
forsterite. In the hydrodynamics module, monomers and dust densities are advecte d wi th the gas density, with the 
terminal velocities given by the low-Reynolds-number case of Eq. (19) in Rossow IU48I1 as settling speed added to 
the vertical advection velocity of dust grains. In contrast to the sophisticated treatment of the gas opacities, a simple 
formula for the dust opacities is used, which assumes that the large-particle limit is valid for all grain sizes and treats 
scattering as true absorption. The dust opacity in each cell of the simulated atmosphere is added to the gas opacity 
(Sect. l3.6T2l . Experiments have been made with another dust model, that uses in addition to one density field for the 
monomers a number of further fields, one for each possible grain size. 



4. Results 

4.1. Code comparison: the solar benchmark 

The natural benchmark for the comparison of different codes is of course the solar atmosphere. On the one 
hand, its mean thermal stratification is well known empirically, and its velocity field and associated temperature 
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Table 1: Setup and emergent radiation of solar models computed with different codes. 





STAGGER 


MURaM 


C05BOLD 


C05BOLD 








standard 


high resolution 


Box size [Mm 3 ] (x X y X z) 


6.0x6.0x3.6 


9.0x9.0x3.0 


5.6x5.6x2.3 


5.6x5.6x2.3 


Grid dimension 


240x240x230 


512x512x300 


140x140x150 


400x400x300 


Cell size [km 3 ] (Ax X Ay X Az) 


25.1x25.1xAz* 


17.6x17.6x10.0 


40.0x40.0x15.1 


14.0x14.0x7.5 


Height range [Mm] (z = at <r> = 1) 


-2.72... + 0.88 


-2.00... + 1.00 


-1.38. .. + 0.88 


-1.38... + 0.88 


upper boundary condition 


transmitting 


closed 


transmitting 


transmitting 


lower boundary condition 


open 


open 


open 


open 


# snapshots used 


19 


19 


19 


60 


# frequency groups used 


12 


4 


12 


12 


Effective temperature [K] 


5762 


5768 


5781 


5763 


Bol. intensity contrast (/i-1) [%] 


14.9 


15.4 


14.4 


14.1 



* The STAGGER code uses a non-equidistant grid in the vertical direction, with spacings ranging from Az — 1 km 
near the optical surface to Az = 32 km in the deepest layers. 



fluctuations have been studied in great detail based on a large body of observations. On the other hand, many numerical 
simulations have been carried out to study solar surface convection with a variety of different computer codes. Here, 
we compare some basic quantities obtained from numerical simulations of the solar surface layers with three different 
codes: STAGGER, MURaM, and C05BOLD. All three codes solve the time-dependent equations of compressible 
(magneto)hydrodynamics for a gravitationally stratified, radiating fluid in a Cartesian box in 3 spatial dimensions, 
taking into account partial ionization and non-gray radiative energy exchange, the latter treated with the opacity- 
binning scheme (see Sect. l3.6.2l ). 

The codes have been developed independently and use different numerical methods. STAGGER and MURaM 
are similar in that both use a method of lines for the hydrodynamics part as well as artificial (hyper)diffusivities to 
stabilize the numerical solution. The STAGGER code 115911 . 114911 (see also Sect.[TJ uses a sixth-order finite-difference 
method to determine the spatial derivatives on a staggered mesh, while the equations are stepped forward in time using 
an explicit third-order predictor-corrector procedure, conserving mass, momentum, energy, and magnetic-field diver- 
gence. Ra diativ e energy exchange is found by the formal solution of the Feautrier equations on long rays. Similarly, 
MURaM lfl5(ill . 15411 uses a fourth-order central-difference scheme in space, and a fourth-order Runge-Kutta time 
stepping; radiation transport is computed with a short-characteristics method. On the other hand, C05BOLD is based 
on a finite-volume approach and employs an approximate Riemann solver of Roe type to advance the hydrodynamics 
in time, relying on second-order monotonic reconstruction schemes to achieve numerical stability without the need 
to invoke artificial viscosities. Directional splitting reduces the 3D problem to one dimensional sub-steps. Similar to 
STAGGER, radiative transfer is treated with a Feautrier method on long characteristics (see Sect. l3.6.3l . 

The basic setup of the different solar simulation runs is summarized in TableQ] The two C05BOLD models differ 
only in their spatial resolution. Since the models of the different groups have not been constructed for the purpose of 
this comparison, they differ in many aspects, such as horizontal box size, vertical extent, spatial resolution, boundary 
conditions, opacity tables, number of frequency groups, and equation of state (EOS), apart from the different numerical 
methods used to solve the equations of hydrodynamics and radiative transfer. Despite these substantial differences, 
the mean vertical structure, obtained from the various simulations by horizontal and temporal averaging, turns out to 
be remarkably similar, as demonstrated in Fig,|9] Obviously, the mean thermal structure is the most robust quantity. 
Ignoring the layers influenced by the top boundary, the temperature differences are everywhere below 2%; deviations 
seen in the deeper layers are probably related to differences in the EOS. Except for the photospheric layers above 
sa 300 km, where the details of the opacity-binning recipe play a major role, the predicted amplitude of the horizontal 
temperature fluctuations is also amazingly similar. As a consequence, the predicted continuum-intensity contrast (see 
Sect. l4.2t is found to be in very good agreement (last row of TableQ}. 

The depth-dependence of the mean vertical velocity obtained from the three different simulations agrees closely 
(lower set of curves in lower right panel of Fig,|9]l- As theoretically expected, (V z ) is positive in the convectively 
unstable layers below the surface, and negative in the overshoot region. Somewhat larger deviations among the 
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Geometrical height [Mm] Geometrical height [Mm] 



Figure 9: Comparison of the average vertical temperature structures, (T)(z) (left panels), the rms horizontal temperature fluctuations ST^s = 

yl(T 2 )(z) - (T) 2 (z) (top right), and the mean and rms vertical velocity (V-)(z), and -\J(V 2 )(z), respectively (lower and upper set of curves in bottom 
right panel, respectively), as obtained with different codes for the solar simulations described in TablefT] STAGGER (dashed), MURaM (dotted), 
C05BOLD standard (solid). Here (.) denotes averaging over horizontal planes of the numerical grid (constant geometric height z) and over selected 
snapshots in time. 



different models are found in the velocity dispersion, J{V})(z) (upper curves). Is seems that both the location of the 
lower boundary and the spatial resolution have some influence on the resulting velocity amplitude. Nevertheless, the 
overall agreement is very satisfactory. 

We have to keep in mind that the different codes are largely based on the same physical assumptions and approx- 
imations. It may therefore not be too surprising that the resulting atmospheric structures are similar. And it does not 
prove that all details of the models are physically correct. 

The role of the spatial resolution is illustrated in Fig.[lO] where we compare two C05BOLD models that differ 
only in spatial resolution: in the high-resolution model (cf. Fig.Q]), the horizontal cell size is reduced by a factor 2 y2 
with respect to the standard C05BOLD model, while the vertical cell size is reduced by a factor of 2. The mean 
thermal structure is practically unchanged, as is the amplitude of the T-fiuctuations up to the mid photosphere. As a 
consequence, the intensity contrast is not significantly affected by the increased grid resolution (see TableQ]). How- 
ever, in the upper photosphere above z ~ 300 km, the amplitude of both the temperature and velocity fluctuations 
increases somewhat with increasing spatial resolution. The question whether the moderate spatial resolution of the 
standard hydrodynamical models is fully sufficient to account for the "turbulent" character of the solar photosphere, 
and hence for correctly capturing the non-thermal Doppler broadening of spectral lines, is currently under investiga- 
tion. Obviously, this is an important issue in the context of accurate chemical abundance determinations based on 3D 
model atmospheres (cf. Sect. l4.3l >. 
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Figure 10: Comparison of the same quantities as in Fig.|9] but for two C05BOLD solar simulations which differ only in spatial resolution: 
C05BOLD standard (solid) and C05BOLD high resolution (dashed). The thin curves in the upper left panel refer to (T)(z) ± 6T rms (z). 



4.2. Granular intensity contrast 

The granulation pattern visible at the solar surface is a manifestation of convection in the sub-photospheric layers: 
bright granules correspond to hot rising gas, while the dark intergranular lanes consist of cooler downward-sinking 
material. The relative continuum-intensity contrast, 



^{I(x,y,t) 2 ) x , y -{I{x,u,t))l y 
Wx,y,t)) x .ij 



(73) 



of this granulation pattern and its variation from the centre of the solar disk to its limb are important tests for the degree 
of realism of numerical models. For many years, the values deriv ed from observations were significantly lower than 
those calculated on the basis of num erical simulations [e.g. 



151| 



Recently, it was s hown that synthetic continuum- 
intensity maps based on C05BOLD lll52ll and also the code by Nordlund and Stein lll49ll can i ndee d reproduce the 
empirical values quite well, if the instrumental image degradation is taken into account properly IU53I1 . The necessary 
image reconstruction is a very demanding task, which on the other hand turns out to be cr ucial, as was shown for 
obser vations obtained with the Solar Optical Telescope (SOT) onboard the Hinode satellite 111 5411 . Hirzberger et al. 
1 15511 find very good agreement between the rms contrast of solar granulation obtained from measurements with a 
balloon-borne 1-m solar telescope and simulations at wavelengths of 388 nm and 312nm. At shorter wavelengths, 
discrepancies between observations and simulations seem to persist. 

By analogy, the surface of cool stars must be covered by a similar pattern, the stellar granulation. Its intensity 
contrast cannot be measured directly. Numerical simulations are necessary to infer how the intensity contrast depends 
on the stellar parameters: the effective temperature r e ff, the surface gravity log g, and the chemical composition. 
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Figure 1 1 : Top: grid of C05BOLD models in the log T c ff-log g diagram. Bottom: logarithm of rms bolometric intensity contrast versus log T e g. 
The squares mark 3D local models, the triangles 2D local models. The stars indicate global models (with very low surface gravity). Larger symbols 
indicate lower gravity. The Sun has its own standard symbol O with yellow background. A red symbol (at lower effective temperatures) shows that 
some treatment of dust is included in the simulation. 



The top panel in Fig.Qj] shows the C05BOLD models on the main sequence and above it (i.e., with gravities 
around log g 4 to 5 and lower: white dwarf models are not plotted) in a log T e g - log g diagram. The displa yed m odels 
comprise the solar-metallicity part of the CIFIST grid of solar-lik e 3D models 1 15611 . 3D M-dwarf mode ls 1 15711 . local 
2D "dusty" brown dwarf models l96ll . global 3D red supergiant 1 158 , 159 , 160ll and AGB star models lll44ll . as well 
as more experimental models of e.g., A-type stars (in 2D or 3D) and cepheids (in 2D). Larger symbols mean lower 
gravity (and usually a larger stellar radius). Squares depict 3D models, triangles 2D models. Solar models have the 
O symbol. Global 3D models of red supergiants and AGB stars are marked as star symbols at the top. Red symbols 
indicate that the simulations have accounted for dust in some form. Models with non-solar metallicities are not shown. 

The bottom panel shows the (bolometric) relative intensity contrast according to Eq. (1731 versus T e ff for the same 
models and with the same symbols as in the top panel. On the main sequence (smallest symbols), the contrast 
decreases for stars cooler than the Sun since the stellar energy flux decreases and convection can transport it with 
smaller temperature fluctuations. The contrast does not increase further but has a plateau for stars a bit hotter than 
the Sun because convection does not carry the entire stellar flux anymore. Below a minimum at around 2600 K, the 
contrast increases again because fluctuating dust clouds start dominating the surface contrast (see Sect. l4.6l l. The 
contrast decreases at the very cool end due to the decreasing overall flux. 
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In general, lowering the gravity has a similar effect as increasing the effective temperature, but results in slightly 
more vigorous convective flows. The largest surface contrast is seen in the global AGB-star models, followed by the 
global red-supergiant models. 2D models have a larger contrast than 3D models. Other types of dust (and/or dust 
schemes) as well as global fluctuations might change the picture for the cooler models. 



4.3. Solar and stellar abundances 



One important application of multi-dimensional (multi-D) radiation-(magneto)hydrodynamics models is the de- 
termination of chemical abundances in late-type stars. Under most circumstances, information about the thermal 
and kinematic structure of a stellar atmosphere is necessary to interpret measured strengths of spectral lines in terms 
of chemical abundances. To this end, simulated time series of the evolution of the stellar photospheric flow field 
are serving as input for detailed spectral-synthesis calculations. The result of these calculations are time series of 
spatially resolved synthetic spectra which, after suitable averaging in space a nd tim e, can be compared to the obser- 
vations. For C05BOLD, we developed the spectral-synthesis code Linfor3D Il6lll which is on the one side adapted 
to the particular data formats and structures of C05BOLD, and on the other side designed to facilitate the abundance 
analysis. 

Historically, the application of multi-D models for deriving abundances started out on the Sun already early on. 
However, in the beginning mainly structural properties of surface convection and associated magnetic fields were in 
the modeling focus so that abundan ce analyse s with multi-D models remained sparse. A turning point came with the 
work of Asplund and collaborators Il62lll63ll suggesting that multi-D effects are important in the Sun and metal-poor 
stars if one wishes to obtain high-fidelity abundances. Since then efforts are directed towards improving multi-D 
modeling aspects specific to abundance ana lysis work, and extending the model basis covering successively larger 
regions of the Hertzsprung-Russell diagram 115611 . 

Hitherto, C05BOLD models were applied to derive abundances of twelve elements in the Sun (see Caffau et al. 
f33ll and references therein) including the CNO elements, which are important for the overall solar metallicity; work 
on further elements is in progr ess. In the field of metal-poor stars, CQ5BOLD models were used to derive abundances 
from atomic (e.g., 1 164 . 165 1) as well as molecular lines (e.g., 1 166l 167]). An element of particular interest in 
metal-poor stars is lithium due to its connection to nucleosynthetic processes in the big bang and early universe. The 
lithium abundance is commonly derived from the Li I resonance line at 6707 A. Since lithium is mostly ionized in 
the stars of interest, the formation of the line is highly temperature sensitive, which makes the resulting abundances 
strongly model-dependent. Hoping for lithium abundances of higher fidelity, multi-D models were rather exten sivel y 
applied. C05BOLD models were used to obtain lithium abundances in the most metal- poor dwarf stars known II681I . 
to investigate the structure of the so-calle d "Sp ite plateau" at lowest metallicities Il69ll . and to study the evolution of 
lithium in the globular cluster NGC 6397 II 7011 . A related aspect is the abundance ratio between the lithium isotopes 
6 Li/ 7 Li i n metal-poor stars. C05BOLD models were applied to argue against claims of a non-zero isotopic ratio 
I 34lll7lll . The aforementioned investigations focused on dwarf or subgiant stars. There are also o ngoing efforts t o 
extend the application of C05BOLD simulations to giant stars including studies of their abundances Il72ill73lll74ll . 

Besides conducting actual abundance analyses, C05BOLD models were instrumental in a number of studies more 
indirectly linked to the derivation of stellar abundances from spectroscopy: the long-lasting issue of how sm all-scale 
velocity fields in stellar atmospheres give rise to the spectroscopically derived microturbulence 1 175l 157 1. and the 
influence of thermal inhomogeneities on effective temperatures derived from Balmer-line profiles lll76ll . 

The application of multi-D models to stellar-abundance studies is in its early stages and modeling challenges still 
exist: the need for a precise thermal structure of the optically thin, line-forming regions demands for a detailed rep- 
resentation of the radiation field. Sufficient wavelength resolution (Sect. [3~l6~2"l i. and inclusion of scattering processes 
are current challenges in the simulations proper. In the post-processing, line-formation calculations including depar- 
tures from LTE are necessary to fully exploit the model potentialities but are demanding in terms of computational 
resources and amount of necessary atomic input data. 



4.4. The magnetic Sun 
4.4.1. Current status 

Figure[T2]exemplifies the type of MHD applications that are currently performed with C05BOLD. It illustrates the 
magnetic field structure at the interface between the convection zone and the overlying atmosphere of the Sun. The top 
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Figure 12: Flux expulsion in a close-up from a simulation of solar magnetoconvection: Logarithmic magnetic field strength in a vertical cross 
section (top) and in three horizontal cross sections (bottom) at heights of km, 250 km, and 500 km. The emergent intensity is displayed in the 
rightmost panel. The arrows represent the velocity field projected in the respective coordinate plane s. Th e yellow/white curve in the top panel 
marks the height of visible optical depth unity, i.e., the "solar surface". From Wedemeyer-Bohm et al. Il77ll . 

panel shows a close-up of a vertical cross section through the three-dimensional computational domain, where colors 
represent the magnetic field strength and arrows the velocity field. The dashed yellow/white curve indicates optical 
depth unity, i.e., the "solar surface" as seen in the visible part of the spectrum. Below this surface, the atmosphere is 
convectively unstable and energy is transported mainly by convection. Above this surface adjoins the stably stratified 
photosphere, where energy is mainly transported by radiation and where convective overshoot motions are rapidly 
damped. The panels in the bottom row show horizontal cross sections of corresponding size at three selected height 
levels and the emerging bolometric intensity in the rightmost panel (intensity map). The location of the vertical cross 
section is indicated by the dashed horizontal line in the bottom panels. 

We can see a strong central updraft in the vertical cross section of Fig.[T2] which corresponds to the central granule 
visible in the intensity map. This granule is a typical representative for real solar granulation with respect to intensity 
contrast and size. Since the diffusion length scale of the magnetic field is small compared to the size of a granule, it is 
useful to think of the magnetic field to be "frozen into the plasma" so that the flow field transports the magnetic field 
from the granule center to its boundaries where it gets concentrated. This process is called the flux expulsion process, 
as magnetic flux is expelled from the granule interior to its boundaries. Correspondingly, the magnetic field in the 
central part of the granule is weak (dark blue) and it gets concentrated in the intergranular lanes (red), where plasma 
flows back into the convection zone again. 

As the updraft runs into the stable stratification of the photospheric layer, it loses the driving buoyancy force and 
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Figure 13: Logarithmic current density, log in a vertical cross section (top panel) and in four horizontal cross sections (bottom panels) in a depth 
of 1 1 80 km below, and at heights of 90 km, 610 km, and 1310 km above the mean surface of optical depth unity from left to right, respectively. The 
arrows in the top panel indicate the magnetic field strength and directio n. Th e dashed line in the bottom row indicates the position of the vertical 
section, j is given in units of 3 X 10 5 A/m 2 . From Schaffenberger et al. Il83ll . 



buoyancy starts to act in the opposite direction. Consequently, and also because of the strong density stratification, the 
flow must deviate in the horizontal direction and it drags the magnetic field with it as a consequence of the frozen-in 
condition. Hence, the magnetic field assumes a predominantly horizontal direction in the upper part of the photo- 
sphere, above the mushroom shaped void in Fig. [12] (top panel). While MHD simulation have since long predicted 
the existence of this prevalently horizontal field 117811 . it is only v ery r ecently that it got observationally detected 
with polarimetric measurements from the Hinode space observatory J179D. More details abo ut MHD simulation s wit h 
regard to this horizontal field can be found in Schiissler and Vogler 118011 . Steiner et al. 118111 . and Steiner et al. 1118211 . 

FigureQj] shows the electric current density, j = V x B, that forms as a consequence of the process discussed 
above, in a similar but larger domain than that of Fig. [12] Again, the top panel is a vertical cross section through 
the simulation domain, and the bottom row shows four horizontal cross sections at various heights. The magnetic 
flux concentrations in the intergranular lanes, which form as a consequence of the flux expulsion process, give rise 
to the conspicuous pairs of current sheets, visible in the horizontal cross section at z — 90 km. Another system of 
current sheets forms in the region of predominantly horizontal magnetic fields above granules, in a height range from 
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approximately 400-900 km, as can be seen from the vertical cross section. Higher up in the atmosphere, shock waves 
form where the supersonic plasma flow sweeps magnetic field into the compression zone downstream of and along 
the shock front. There again, current sheets form, e.g., the one extending from x = 1.1 to x = 1.5 Mm in the top 
part of the vertical cross section, which is also visible as a thin filament at x = 1.4 Mm in the rightmost panel which 
corresponds to z = 1310 km. More details about MHD simulations with CQ5B OLD with regard to chromospheric 



shock waves can be found in Schaffenberger et al. 1530 and Schaffenberger et al. 1183 1 

Interestingly, the mean vertical Poynting flux, (S 3), where S = B x (v x B) = B 2 v - (v • B)B is all of the magnetic 
contribution to the total energy flux (see Eq.(fT9ll), changes sign near optical depth unity. Below this depth, in the 
convection zone, the intense downdrafts in the intergranular lanes pump magnetic fields in the downward direction — 
cool and dense plumes compress and drag the magnetic field with them. This leads to a net Poynting flux in the 
downward direction. Above t c = 1, the Poynting flux that is connected to the magnetic field carried by the convective 
overshoot prevails, leading to a net Poynting flux in the upward direction. We expect that at least a part of this flux 
is turned into heat via ohmic dissipation by the current sheets that form on top of the overshoot in the chromospheric 
layers (see Fig. [TBI. 

One purpose of performing realistic simulations is the synthesis of observable quantities from the simulation data, 
which then become directly comparable to actual observations of the Sun. This typically involves the integration of the 
radiation transfer equation along lines of sight across the computational domain in order to obtain (two-dimensional) 
synthetic intensity maps. This analysis step is performed post factum, after completion of the simulation. In case 
of magnetohydrodynamics, it requires the integration of the Unno-Rachkovsky equation for polarized radiation for 
obtaining intensity maps of the Stokes parameters. For direct comparison with observations from space-based or 
ground-based observatories, application of the corresponding instrumental point spread-function is necessary for de- 
grading the synthetic intensity maps to the spatial resolution limits of the actual observation. Degradation in frequency 
space and addition of noise may be required for taking into account the frequency resolution of the spectrometer and 
the photon noise of the recording device, respectively. Quantitative comparisons of synthesized spectropolarimetric 
maps from C05BOLD sim ulatio n data with corresponding intensity maps from the Hinode space observatory were 



perfo rmed by Steiner et al. 118111 who focus on the above discussed horizontal magnetic fields and by Rezaei et al. 
1 18411 with respect to intergranular magnetic flux concentrations. 

MHD simulations with C05BOLD are also performed for studying the excitation and propagation of magneto- 
acoustic waves in the magnetically structured solar atmosphere. Effects of mode conv ersion, refraction, and transmis- 
sion are studied for application in solar atmospheric seismology 11851 . fl86lfl87lfl88ll . 



4.4.2. Next steps 

Most if not all simulations of stellar magnetoconvection rely on the MHD approximation. But the plasma in the 
photosphere of the Sun is weakly ionized so that the frozen-in condition ma y not real ly apply after all, despite the 



large scales of the magnetoconvective processes under consideration [see, e.g., fl89l[l90ll ." In the tenuous plasma of the 
chromosphere and corona, effects become important which are not included in the single-fluid model of the standard 
MHD equations, and a multi-fluid model or even a kinetic description of the plasma may become necessary. 

Two consequences of a multi-fluid description are the Hall effect and the ambipolar diffusion. Since the electrons 
and ions are moving on curved trajectories between collisions, the current density vector is no longer colinear with the 
electric field vector, i.e., the electric conductivity is a tensor. Ambipolar diffusion occurs in partially ionized plasmas, 
where collisions between charged and neutral particles produce new diffusion mechanisms. Since only the charged 
particles are coupled to the magnetic field, the forces acting on the different particles are different. This leads to a drift 
between charged and neutral particles modifying the transport of magnetic flux in the plasma. 

In the limit of a weakly ionized plasma consisting of electrons, ions, and neutrals, the single fluid description 
can still be applied along with an appropriate modification of the induction equation. This modification includes a 
term representing the Hall effect and a term representing the ambipolar diffusion. We plan to implement an optional 
inclusion of these terms in the computation of the numerical fluxes for the induction equation of C05BOLD. 

Ohmic dissipation of electric currents is probably the most important process in the heating of the outer solar 
atmosphere. In the present version of the MHD module of C05BOLD, there is only an explicit implementation of a 
turbulent subgrid-scale magnetic diffusion (see Sect. 13.7.51 . However, for taking a significant resistivity into account, 
an implicit treatment must be implemented. Also the occurrence of anomalous resistivity should be accounted for. 
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Currently, a much discussed topic in solar physics is the origin of the ubiquitous weak magnetic field outside 
sunspots. It is thought to be generate d thro ugh induction due to the turbulent motion of the plasma near the solar 
surface, viz., by the turbulent dynamo [191]. On the other hand, turbulent pumping transports this magnetic field in 
the downward direction, away from the surface into the deep, less turbulent layers of the convection zone [192]. It is 
unclear, which effect prevails. First MHD simulations have been designed and carried out with the aim to improve 
our understanding of this riddle (Vogler and Schiissler B193I1 . Moll et al. I194l0 . For a review see Stein 1 195] who 
argues that the solar dynamo has no preferred scale but rather acts throughout the convection zone over a wide range 
of scales. 

Regarding global models of solar magnetohydrodynamics, the long term goal should be the global, three-dimensional, 
numerical simulation of the entire solar convection zone including at least the overshoot layers and a proper radia- 
tive transfer at its boundaries, but ideally also photosphere, chromosphere (see below), and corona. Thus, we seek a 
"Sun simulator" as a laboratory for the holistic simulation of the Sun on scales from the stellar radius to the size of 
granules of about 1000 km. The outcome should be a virtually self-consistent simulation of the differentially rotating 
convection zone including the self-exciting dynamo, which is thought to be at the origin of solar magnetic activity. 
Simulations should shed light on the functioning of the dynamo, the solar rotation law, the torsional oscillation, the 
luminosity variability, the sunspot cycle, and the global solar oscillation. Despite the apparent sphericity of stars, these 
processes are truly three-dimensional and therefore, they require a three-dimensional treatment unlike the traditional 
one-dimensional approach in stellar-evolution modeling. 

The development of such a simulation tool is a formidable task (see Sect l2.U . The main challenge at the beginning 
consists in the recognition of the relevant physics to be included and in finding the proper physical approximations, 
numerical scheme, and adaptive meshing for achieving sufficient spatial resolution. So far, global MHD simulations 
have not been carried out with C05BOLD. However, with C05BOLD it should be possible to collect first experiences 
on the way to a true "Sun simulator". 



4.5. Solar chromosphere 

The chromosphere is the thin atmospheric layer between the photosphere and the transition region and corona 
above. Although these layers are coupled by magnetic fields and waves, the properties of the atmospheric gas differ 
significantly. Compared to the photosphere below, the chromospheric gas is relatively thin, which has a number of 
important implications for the modeling (see Sect. l4.5.TI ) and also for observations. There are only a few diagnostics 
suitable for probing the chromosphere, which makes it hard to derive constraints and viable reality checks for numeri- 
cal models. Even worse, the interpretation of most of these diagnostics is complicated by the fact that non-equilibrium 
effects must be taken into account. Advances in instrumentation both for ground-based and space-borne observations 
during the recent years made it nevertheless possible to access the dynamic and intermittent fine structure at small 
spati al scales like it is seen in current chromosphere simulations (see, e.g., the review by Wedemeyer-Bohm et al. 
117711 ). The coexistence of magnetic fields and propagating waves and interaction of these makes the modeling of the 
chromosphere a challenging task and a true hardness test for the stability of the code. 



4.5.1. Challenges in chromospheric modeling 

Radiative transfer: The gas becomes optically thin (i.e., essentially transparent) in a substantial wavelength range, 
leading to a strongly non-local coupling of regions within this layer and also with the layers below and above. How- 
ever, the chromosphere is neither completely optically thin nor completely optically thick. The often used simplifying 
assumption of local thermodynamic equilibrium (LTE) breaks down in the chromosphere and effects like scattering 
become important. All this makes the numerical treatment of radiation challenging. Detailed (non-LTE) calculations 
are today usually done for (a number of) single simulation snapshots but are still computationally too expensive to be 
included in 3D radiation magnetohydrodynamic simulations. Simplifications are unavoidable, so far. 

Non-equilibrium effects: Many more processes depart from equilibrium conditions in the thin chromosphere. For 
instance, the ionization degree of hydrogen can no longer be modeled under the assumption of an instantaneous 
equilibrium (see Sect. 13.8.2b . Some processes become so slow that the detailed ti me evolution m ust be followed. 
A detailed treatment, however, is computationally expensive. Current approaches 1 140i 141 , 196 1 use a hydrogen 
model atom with 6 energy levels and 10 radiative transitions. The corresponding rate matrix can be stiff and is 
solved implicitly for each grid cell in the model chromosphere. Another example for non-equilibrium modeling is the 
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application of chemical-reaction networks for carbon monoxide 11 351, 1 1371 1 19711 . See Sect. l3.8l for a short description 
of the implementation of chemical-reaction networks and hydrogen ionization in C05BOLD. 

Computational time step: The thermal pressure is so low that the magnetic pressure can become larger already at 
relatively low magnetic field strengths. The plasma-/?, which is defined as the ratio of thermal to magnetic pressure, is 
consequently less than one above heights of ~ 1000 km above the bottom of the photosphere. Under these conditions, 
the magnetic field is no longer advected passively and magnetic wave modes become important. A full magneto- 
hydrodynamic approach is therefore necessary. The computational time step is then determined by the Alfven speed, 
which easily can result in steps of the order of milliseconds (in simulations that cover a few hours), depending on 
the magnetic configuration in the chromosphere. This is a reduction by one to two orders of magnitude compared 
to purely hydrodynamic simulations. Under certain circumstances, an artificial reduction of the Alfven speed can be 
used, as is discussed in Sect. 13.7.41 

Numerical stability: Acoustic waves, which are continuously excited by the non-stationary surface convection 
below, grow in amplitude while propagating into the thinner chromosphere. There, they develop (MHD) shocks with 
high peak temperatures of the order of 7000 K or more, and the dynamical pressure exceeds the gas pressure. The 
physical conditions in a grid cell can change drastically during the passage of a shock wave, which requires a high 
degree of stability of the numerical scheme. In the MHD case, the occurrence of strong gradients in thermal and 
magnetic properties can lead to exceptional situations in small parts of the computational domain. In this respect, the 
HLL solver (see Sect. 13. 7b has been proven a good choice. 



4.5.2. Chromospheric modeling in the recent years 

The many complications in modeling the chromosphere made it inevitable to begin with simplified models and 
increase th e degree of realism step by step. A very prominent example is the pioneering study by Carlsson and Stein 
They restricted the sim ulatio ns to one spatial dimension but implemented a detailed radiative transfer 



treatment [cf. 120011 . Skartlien et al. 1120 111 succeeded to produce a 3D hydrodynamic model with a relatively coarse 
spatial resolution by using a simplified description of the radiative transfer, which nevertheless included scattering. 
The 3D hydrodynamic simulations by Wedemeyer et al. 119211 . which were carried out with C05BOLD, did not include 
scattering but used a higher spatial resolution. This type of local 3D models is restricted to a relatively small part of 
the atmospheric layers extending from the chromosphere into the upper convection zone. The latter is important as 
it provides an intrinsic driver for the atmospheric dynamics and thus avoids the need for an artificial driver like it is 
necessary in ID simulations. 

The step to multi-dimensional magnetohydrodynamic simulations of the chromosphere has been performed only 
a few years ago. The (non-local) radiative transfer is still a limiting factor. Consequently, simplifications of the 



radiative transfer are still necessary for 3D MHD simulations of the solar chromosphere. Schaffenberger et al. 15311 



therefore used a frequency-independent ("gray") radiative transport and a weak initial magnetic field for the fir st 3D 



MHD simulations with C05BOLD (see also Schaffenberger et al. [183])- Further 2D numerical experiments [202] 
combined higher magnetic field strengths (Bo = 100G) with the treatment of chemical -reaction networks including 
carbon monoxide and the methylidyne radical CH in view of their diagnostic potential. 

The models mentioned above focus on the small-s cale struct ure and dyn amics of the solar chromosphere, while 



another class of models also includes the corona above 1203L 1204 l65l l64l I205L 15211 . These models have a larger spatial 



extent so that the large-scale magnetic field structure can be fitted into the computational box. In order to keep the 
simulations feasible, compromises such as a lower spatial resolution were unavoidable for the earlier models. How- 
ever, the progress in computational performance and efficient numerical methods allows for higher spatial resolution 
and at the same time a larger number of implemented physical processes, producing models with a increasing degree 
of realism. 



The chromospheric layer of the hydrodynamic C05BOLD models by Wedemeyer et al. 119211 exhibits a very dy- 
namic and intermittent pattern made of propagating hot shock fronts and cool post-shock regions (cf. Fig.[T4"b). The 
resulting fluctuations of the gas properties are substantial like it was found already from ID simulations. As the shock 
fronts are very narrow, the peak temperatures of 7000-8000 K in C05BOLD simulations depend to some degree on 
the resolution of the computational grid. Adiabatic expansion of the post-shock regions produces gas temperatures 
down to ~ 2000 K. A similar shock pattern can already be perceived in the gas-temperature maps by Skartlien et al. 



1 201] but much less clearly due to the lower spatial resolution. Martmez-Sykora et al. [64], who employed the Oslo 
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Figure 14: Vertical cross sections through the upper layers of a 3D MHD model by Wedemeyer-Bohm .206 1. a) logarithmic magnetic field strength 
(color range 0.2 - 3.5 [G]), b) gas temperature (color range 2500 - 12000 K). The solid curves in both panels represent (projected) magnetic field 
lines (thin black curves), the contour for plasma-/? = 1 (thick yellow contour), and the height where the optical depth is unity (horizontally running 
black curve around z = 0). The latter defines the bottom of the photosphere. 



Staggered Code for 3D simulations of magnetic-flux emergence, find a shock-induced pattern in their model chro- 
mosphere, too. The range in chromospheric gas temperature is similar to the C05BOLD results, whereas there are 
differences in the temperature amplitudes. This is presumably caused by the different numerical treatment of the ra- 
diative transfer in the upper layers. The exist ence of a c hromospheric small-scale pattern in quiet regions of the Sun 
is now supported by recent observations [e.g. 



2071 2081. 



Also the MHD simulations carried out with C05BOLD exhibit strong shock fronts in the chromosphere (see 
Fig. [14b). Compared to the photosphere below, the magnetic field in the model chromospheres is le ss concentrated 
and reaches a higher filling factor; it has a lower average field strength and evolves faster II53U 18311 . See Sect. 14.41 



for a description of the photosphere in this type of models. The topology of the chromospheric magnetic field is yet 
complex and features shock-induced compression and amplification into magnetic field filaments. Fig.[T4"a. shows the 
upper layers of a current 3D simulation. The vertical cross section is intersecting a magnetic flux concentration. The 
magnetic field lines (thin solid lines, projected into the view plane) show that the field is highly concentrated in the 
photosphere and expands in the chromosphere above. The thick yellow curve represents the surface where plasma- 
P = 1. Although the height of this surface varies strongly, it is typically found around z ~ 1000 km outside strong 
magnetic flux concentrations. For plasma-/? = 1, sound speed and Alfven speed are sim ilar, which has importan t 
implications for the occurrence, conversion, and propagation of different wave modes [e.g., l27^ElHEnlfl85[fl87M . 
The simulations indeed show a different behavior in the domains with j3 < 1 and j3 > 1, resulting in a slowly evolving 
lower part and a highly dynamic upper part. Current sheets exist below and above the [3-1 surface but differ in 
their orientation (see Fig. [TBI. They are stacked with predominantly horizontal orientation in the lower atmosphere 
(see Sect. l4.4b but are aligned with shock fronts in the low-/? regime in the chromosphere, resulting in oblique or even 
vertical orientation. 



4.5.3. Next steps 

The next steps towards realistic models of the solar chromosphere concern an improvement of the radiative transfer 
under chromospheric conditions and the detailed treatment of non-equilibrium processes that have a significant impact 
on the equation of state and the opacities. However, the detailed modeling of non-equilibrium effects might increase 
the computational costs to a degree that it can become impractical. For instance, important opacity sources that deviate 
from their equilibrium state can prevent the usage of numerically efficient opacity look-up tables (Sect. 13.6.21 ) and 
require a costly detailed line-by-line treatment of the radiation transfer. The inclusion of large-scale magnetic fields is 
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another important point as most of the models discussed above resemble rather quiet Sun internetwork conditions with 
comparatively weak magnetic fields and thus cannot be applied to more active regions on the Sun. Chromospheric 
simulations with C05BOLD have been restricted to heights below the transition region, where thermal conduction 
can still be neglected. Simulations that include the transition region and low corona have not been performed so far 
because it would require the implementation of the computationally expensive treatment of thermal conduction. 



4.6. Local models of other stars 

Surface-intensity snapshots from the hot A-type star regime (ab out 8500 K) ove r solar-type stars (r e ff s U n=5775 K) 
to brown dwarfs (about 1500 K) are shown in Fig. [15] A-type stars 112 12ll213ll2 1411 have only thin surface convection 
zones, where convection carries only a small fraction of the total energy flux. Of interest in these stars is the amount 
of overshoot below the convection zone, that causes a mixing of elements counteracting the separating eff ect o f 
gravitational settling and radiative acceleration l36ll . The shapes of spectral lines differ from solar counterparts 1 214 1. 
suggesting different atmospheric flow patterns or deviating correlations between temperature fluctuations and velocity 
fields. Due to the shorter radiative time scales, caused by efficient radiative energy exchange, the steeper and stronger 
subphotospheric temperature jump, and the larger convective cells, the simulations are numerically more challenging 
than those of solar-type stars, requiring more numerical grid cells and many more time steps. Therefore, the required 
CPU time per model goes up by a factor of 100 or more - depen ding on the stellar parameters - and an implicit 
treatment (at least of the radiation transport) seems appropriate 121511 . The transition from a thin, inefficient convection 
zone to a deeper zone, where the convection carries in some layers almost all the energy flux, occurs in a similar way 
in the temperature sequence of 2D cepheid models at log g-2 in Fig. [16] Remarkable is the extended overshoot region 
with significant negative convective flux. 

F, G, and K dwarfs form a temperature sequence, in which the convection zone gets deeper, the Mach number 
of the convective velocity declines, the relative efficiency of convection compared to radiation increases, and the 
granular contrast decreases (Fig.[TTl). While the amplitude of pressure waves drops rapidly, the amplitude of gravity 
waves decreases more slowly: in solar models they are more difficult to detect than pressure waves l39ll . but they 
dominate in brown dwarfs 19611 and influence the shape of the dust clouds (coolest models in Fig. [15]). The change in 
the amount of photospheric overshoot and the optical depth at the top of the convection zone affect the appearance 
of granules i n Fig . [131 The scale of the granules seems related to the surface - or rather the sub-surface - pressure 
scale height 121611 . Both are indicated by the horizontal bars in Fig. [15] which have lengths of 10 H p , measured at two 
different heights. 

The first multi-D rad iation-hydrodynamical model atmospheres for M-ty pe st ars were calculated by Ludwig and 
collaborators 1 217 , 218 1 using the simulation code of Nordlund and Stein l219ll . An important issue, that needed 
to be settled in the model construction, was the handling of molecular opacities in the opacity-binning scheme (see 
Sect. 13-6 -"St . It turned out that no particular treat ment is necessary, as long as the molecular opacities dominate the 
atmospheric opacities. Subsequently, Wende et al. IU57I1 used C05BOLD with the previously developed opacity set-up 
to calculate a sequence of models covering the main-sequence in the temperature range 2600 K < T e g < 4000 K, and 
probing surface gravities 3.0 < \og(g) < 5.0 at fixed effective temperature of * 3300 K. Motivated by observational 
demands, the authors investigated the impact of the velocity field and thermal structure on properties of FeH lines. 

While the 2616 K model of an M-dwarf in Fig.[l5jdoes not show visible amounts of dust, in brown dwarfs at even 
lower temperatures, dust clouds begin to form in the atmospheres ll96Tl . that show up as small dark patches (shaped 
by gravity waves) above the low-contrast granules in the 2249 K model in Fig. [15] The coolest model in Fig.fTBIhas 
completely opaque dust clouds, that hide the underlying gas convection zone. A long-wavelength and high-amplitude 
gravity wave creates the large-scale pattern in the plot, while convection within the dust clouds causes the thin bright 
filaments. The numerical challenges in these simulations come from the - hardly known - non-equilibrium dust 
chemistry, the long relaxation times of dust settling and mixing, and the need to account for the likely interaction of 
small "cloud" scales, covered by the local models in Freytag et al. l96Tl . and global "weather front" scales, far beyond 
the size of the computational domain of the local models. 

Cool hydrogen-rich white dwarfs are old stars, that have used up their nuclear fuel and cool down slowly. The 
large surface gravity causes high atmospheric pressures, small convective scales, and a settling of elements heavier 
than hydrogen fr om th e atmosphere into subsurface layers. Recently, first models for those objects were computed 
with C05BOLD 822011 . The work follows up on earlier investigations by some of the authors (HGL, BF, MS) and 
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Figure 15: Frequency-integrated-intensity snapshots of global models of a red supergiant and an AGB star with low surface gravity (top row) and 
local models of stars near the main sequence with larger gravities. The title lines show the effective temperature, the decadic logarithm of the 
surface gravity in cm s 2 , and the relative gray intensity contrast - averaged over a representative time span. The length of the upper bar in the top 
right of each frame is 10 times the surface pressure scale height. The bar below is 10 times the pressure scale height but measured 3 pressure scale 
heights below the other level. 



collaborators 1 109l 221 , 222 1 using a code which was a precursor of C05BOLD. Similar to the case of M-type models, 
the intention is to study the influence of multi-D effects on the formation of Balmer lines. 



4.7. Global models of supergiants and AGB stars 

Early e xplan ations of the irregul ar lig ht curves of red supergiants as due to giant convection cells by Stothers 
and Leung 122311 and Schwarzschild 122411 got more recently support by interferometric detections of spatial inhomo- 
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Figure 16: Convective fluxes (normalized by the flux emitted at the top of the model) versus height for a temperature sequence of cepheid models 
at log g=2. The curves are shifted in z, so that the entropy minimum at the top of the convection zone lies at z=0. 



geneities on the surface of Betelgeuse (e.g., [225]). When convective scales are not very small compared to the stellar 
diameter, global star-in-a-box simulations are possible (Sect. 13. 2^21 making these stars the easiest targets for global 
models, that include a major part of the convective envelope as well as the near environment of the star. The top row 
of Fig.[T5lshows snapshots of the emer gent intensity of such models: on the left from a red-supergiant Il58ll and on 
the right from an AGB-star simulation [ 144]. The supergiant model confirms that convective scales are indeed huge, 
with a few very large, deep, long-lived envelope cells and many small, short-lived surface cells. The surface contrast 
is enormous (see Fig.lTTt. due to violent convective flows and in addition waves, that have already in the lower photo- 
sphere a large amplitude, in contrast to the solar case. Large scales and contrast values render the feat ures observable 
with current interferometers: the models compare favorably with VLTI observations indicating that 

these global models start to become "realistic", too. 

RHD simulations of an AGB star (Fig. [15] top right IU44I0 demonstrate that convection can excite pressure waves 
with amplitudes sufficient to turn them into shocks, which then push out dense material into layers cool enough 
that dust can form (included in the 3D models, see Sect. 13.8.31 ). This allows radiation pressure on dust t o acc elerate 
the material outward causing a stella r wind (no t included yet in the 3D models of Freytag and Hofner 1 14411 but in 
the ID simulations of Hofner et al. 1 106l 107 11. Major challenges for the simulations are posed by the molecular 
opacities varying strongly with frequency, that cause - together with large dust opacities - very small radiativ e tim e 
steps during the time-explicit treatment of the radiative energy exchange. So far, only ID RHD models (e.g., 110711 ) 
include important ingredients like scattering, radiation pressure, and a sufficiently large computational volume to 
account for the extended wind acceleration region. While the properties of the simulated surface granulation seem 
already quite realistic, there are discrepancies further out: The models have a too steep density drop and show no 
"molsphere", chromosphere, or wind. Future generations of models will help to investigate the role of radiation 
pressure on molecular lines and dust, magnetic fields, and rotation for these outer layers. 



5. Conclusions 

For stellar parameters close to the solar values, the transition from ID static stellar-atmosphere models to 3D dy- 
namic local stellar-atmosphere simulations is in full swing. RHD simulations of surface convection of such stars are 
routinely performed by codes like C05BOLD, as presented in this paper. There is a good consistency between the re- 
sults of similar codes and with solar observations. The simulations provide insight into processes related to stellar sur- 
face convection, and deliver high-accuracy model atmospheres for spectrum-synthesis and abundance-determination 
work for a variety of stellar parameters. 

However, many physical effects are not properly incorporated by the current models: small-scale convective struc- 
tures, covered by local-box simulations, interact with their environment via, for instance, large-scale convective flows, 
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magnetic fields, waves, or global dust flows (the latter in cool substellar objects). The inclusion of the chromosphere, 
the corona, the wind-formation zone, etc., requires to cover an even wider range of densities and temperatures than 
the previous models. At lower densities, the detailed treatment of non-equilibrium processes (molecule formation, 
radiation transport not in local thermal equilibrium) are a major challenge, requiring algorithms with a complexity far 
beyond the current treatment. The modeling of magnetic phenomena needs appropriate MHD or plasma-physics sim- 
ulations. Objects significantly cooler than the Sun require a detailed non-equilibrium treatment of dust and "weather 
phenomena". Varying efficiency ratios between radiation and convection as a function of stellar parameters have to 
be considered: the extremely small radiative relaxation time scale in hotter stars causes small numerical time steps 
and slows down simulations significantly. Cool objects on the other hand need extended simulation runs because of 
their long thermal relaxation time scales. Low-gravity objects have extended atmospheres and can produce winds. 
Magnetic-field phenomena exist on very large and very small scales and couple the stellar interior to the photosphere 
and the envelope. In stars more active than the Sun, the fields are stronger and can form very different configurations. 

In the future, we will see a refinement of local simulations and more and more extended model grids, providing 
reliable stellar model atmospheres. However, the main challenges arise from the need to extend the simulations 
in terms of stellar parameters (from A-type stars to planetary objects and from supergiants to white dwarfs), from 
physical effects, and from the extension of spatial and temporal scales towards 3D large-scale or even global dynamic 
models. These models should not only include the photosphere but the stellar interior and the outer atmospheric layers 
as well, covering short and long time scales (many rotation periods, dynamo cycles, stellar oscillation periods, climate 
cycles). 

Realistic global 3D MHD simulations for cool stars will remain a dream for the foreseeable future. Nevertheless, 
numerical simulations will continue to be indispensable tools for the understanding of the various complex dynamical 
processes in stars. 
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